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Abstract 

The Ising model is famous model for magnetic substances in Statistical Physics, and 
has been greatly studied in many forms. It was solved in one-dimension by Ernst Ising 
in 1925 and in two-dimensions without an external magnetic field by Lars Onsager in 
1944. In this thesis we look at the anisotropic Ising model on the Union Jack lattice. 
This lattice is one of the few exactly solvable models which exhibits a re-entrant phase 
transition and so is of great interest. 

Initially we cover the history of the Ising model and some possible applications 
outside the traditional magnetic substances. Background theory will be presented 
before briefly discussing the calculations for the one-dimensional and two-dimensional 
models. After this we will focus on the Union Jack lattice and specifically the work of 
Wu and Lin in their 1987 paper "Ising model on the Union Jack lattice as a free fermion 
model." |WL87j . Next we will develop a mean field prediction for the Union Jack lattice 
after first discussing mean field theory for other lattices. Finally we will present the 
results of numerical simulations. These simulations will be performed using a Monte 
Carlo method, specifically the Metropolis-Hastings algorithm, to simulate a Markov 
chain. Initially we calibrate our simulation program using the triangular lattice, before 
going on to run simulations for Ferromagnetic, Antiferromagnetic and Metamagnetic 
systems on the Union Jack lattice. 
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Chapter 1 

Introduction 



The Ising model is popular in statistical physics as a model for magnetic substances. 
The model can also be applied to systems outside of magnetic substances, such as 
diffusion of a gas and alloying in metals. Initially given to Ernst Ising as his PhD 
project by Wilhem Lenz in 1920, it is still the subject of many studies to this day. While 
Ising looked at the one-dimensional chain |Isi25] . it was extended by Lars Onsager to 
the two-dimensional square lattice model without an external magnetic field in 1944 
|Ons44] . The non-planar two-dimensional model, that is one that is with an external 
magnetic field, and models with higher dimensions, still remain unsolved except in rare 
cases [LCLG94]. 

In this thesis, we will investigate the Ising model on the Union Jack lattice and re- 
entrant phase transitions from the point of view of average magnetisation. In particular, 
we will investigate the results of Wu and Lin from their 1987 |WL87t IWL88] and 1989 
|WL89j papers. For this investigation we use three methods in the thesis. First we 
will perform a theoretical analysis of their results. Next we will develop a mean field 
simulation of the lattice and compare the results to the original exact solution. Finally 
we will perform numerical simulations on the various systems and compare these to 
the theoretical results. 

1.1 Summary of the thesis 

In Chapter [2] we will briefly revise some background information that will be required 
for the rest of the thesis. At the beginning of the chapter the thermodynamic laws 
will be presented and discussed. We will then move on to looking at some basic 
statistical mechanical concepts, such as partition functions and entropy. Once we have 
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a framework from the previous sections we will then move on to a discussion of phase 
transitions. Finally at the end of the chapter we will discuss the topic of magnetic 
substances from a chemistry point of view. 

Using this background theory we will move on to a development of the theory 
of the Ising model in Chapter [3j We will look in detail at the development of the 
model, from magnetic sites, to the one-dimensional model and then extending to the 
two-dimensional square lattice model. At the end of the chapter the results for the 
triangular lattice are given. The triangular lattice is important as it will be used later 
for calibration of the simulation. 

Having developed the theory for these lattices, we will move on in Chapter |4] to look 
at the Union Jack lattice model. Initially in the chapter, we will discuss the work of 
Vaks et al. |VLU66j . They use an older approach to calculate the critical temperatures 
of the system. They also discuss the re-entrant transition, with classification of the 
various phases. Then we will go on to look at the work of Wu and Lin |WL87t rWL89j . 
developing prediction functions for both sublattices. After this we will discuss their 
classification of phases at low temperatures. Once we have finished the presentation, we 
will then perform analysis on interesting systems and discuss the differences between 
the approaches. 

In Chapter [5] we will examine the alternative approach of Mean Field theory. First 
we will review the theory on isotropic lattices giving results for the triangular lattice. 
From this we will develop equations for the Union Jack lattice. These will be in two 
forms, a partially uncoupled set of equations which will be dependent on only one 
sublattice magnetisation and a coupled set of equations. We will present simulations 
results using our mean field predictions and compare them to those from Wu and 
Lin |WL87t IWL89j . We will also discuss the limitations of this approach and any 
disagreements with Wu and Lin's results. The simulation program will be discussed 
later in Appendix [C| 

At this point of the thesis we briefly move away from Ising models, to discuss the re- 
quired theory for the next few chapters. Chapter [6] is concerned with the methods used 
for our numerical simulations. This section will be mainly concerned with Stochas- 
tic processes, in particular the Markov random walk. This will be simulated using a 
Monte Carlo method to reduce the necessary calculations. The Monte Carlo method 
we will be using is the Metropolis-Hasting algorithm. We will adapt this algorithm for 
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the simulation program. The actual program code will not be discussed here, but the 
interested reader may find such a discussion in Appendix [B] 

To calibrate our simulation program, in Chapter[7|we will run simulations on the tri- 
angular lattice. These simulations will range from the simplest isotropic ferromagnetic 
system, to the most general anisotropic antiferromagnetic system. We will compare 
our simulation results to the predictions of Baxter |Bax75j and Stephenson |Ste64j and 
comment on the correlation of results. 

Our numerical simulation results on the Union Jack lattice will then be presented 
in Chapter |8} We will follow the structure of the theoretical analysis in Chapter |4], 
performing simulations on various systems of interest. These simulations will allow us 
to further investigate systems where there was a disagreement between the approaches 
on Vaks et al. |VL066j and Wu and Lin |WL87l IWL89j . Equally we will show how 
our simulation program allows the study of systems that are currently not able to be 
theoretically predicted. 

To conclude the thesis, in Chapter [9] a review of the results from the theoretical, 
mean field theory and numerical approaches will be presented. We will then discuss 
how we can minimise the disagreement effects, such as additional conditions. A list of 
references will then follow. 

1.2 History of the Ising model 

In 1920, Wilhelm Lenz proposed a basic model of ferromagnetic substances to his then 
PhD student Ernst Ising. By 1925, Ising was submitting his dissertation |Isi25j which 
was the first exactly solved the one-dimensional case, and as such later was called the 
Lenz-Ising model. He discovered that there is no phase transition in this case. From 
this result he incorrectly extended it to say that there would be no phase transition 
in higher dimensional cases. Indeed during the early part of the twentieth century it 
was believed by some that the partition function would never show a phase transition, 
as the exponential function is analytic and as such the sum of these functions is also 
analytic. However, the logarithm of the partition function is not analytic near the 
critical temperature in the thermodynamic limit. The thermodynamic limit is reached 
as the number of particles in a system, A^, approaches infinity. So in the infinite system 
the partition function is not analytic. 
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After Lenz and Ising proposed the model, in 1936 Peierls |Pei36j was able to explic- 
itly show that a phase transition occurred in the two-dimensional model. He compared 
the high and low temperature limits and showed that at infinite temperatures all con- 
figurations have equal probability, while at high but not infinite temperatures there 
occurs clumping although the average magnetisation is still zero. At low temperatures 
the model resides in either a state with all the spins being positive or negative. Peierls 
questioned whether it was possible for the system to fiuctuate between these states. In 
his work he managed to establish that the model defines super selection sectors. That 
is, domains that are not connected by finite fiuctuations. 

While greats of physics such as Heisenberg (1928) jHei28j . Kramers and Wannier 
(1941) |KW41al IKW41bj use the model to examine ferromagnetism and properties 
such as the Curie temperature (the critical temperature of substances), it was almost 
20 years before Lars Onsager had presented a solution to the two-dimensional model 
with no external field |Ons44j . His paper and solution, although algebraically complex, 
was a large step in understanding the model further. While other solutions have been 
developed that are less algebraically complex, other than in special cases, it has only 
been solved for a zero external magnetic field by Schultz et al. |SML64j . This solution 
will be presented later in Chapter [3j The derivation of the solution was published 
by Yang in |Yan52j . but Onsager presented the result on the 28 February 1942 at a 
meeting of the New York Academy of Science. 

The extension to three dimensions is more complex. In the paper of Schultz et 
al |SML64j . the method is limited by the transformation of the system into fermion 
annihilation and creation operators. In 1972, Richard Feynman |Fey72| said of the 
three-dimensional Ising model, "the exact solution for three dimensions has not yet 
been found." In 2000, Sorin Istrail |IstOOj . by extending Francisco Barahona's work 
|Bar82] . showed that the solution to the three-dimensional problem as well as the 
two-dimensional non-planar model is NP-complete. He used a method called computa- 
tional intractability, which allows one to discover if a problem can be solved in a feasible 
timefram^ As there are about 6000 such problems, and they are all mathematically 
equivalent, a solution to one would be a solution to all and this is an infeasible result. 
While this result does not remove the possibility of a solution (the P versus NP ques- 
tion is still open), it does mean that the solution is not possible using the methods 
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known today. However, the three-dimensional problem, as with the non-planar two- 
dimensional model, can be solved approximately numerically and in some special cases 

In this section, I have presented a brief history showing the major milestones in 
the development of the model in various dimensions, for the interested reader there are 
many more detailed discussions of the history, such as |Bru67] or |PB94j . 



1.3 Applications of the Ising model 

The Ising model has been applied to a great variety of other physical systems such as 
absorption of gases on to solid surfaces, order-disorder transitions in alloys, concen- 
trated solutions of liquids, the helix-coil transition polypeptides and the absorption of 
oxygen by haemoglobin. |VS95l IMan88[ IBDOBj 

An example, taken from |Man88j . of an order- disorder transition in alloys is the 
binary alloy of Copper and Zinc, brass. Brass has many forms depending on the 
atomic percentage (that is, the number of atoms) of these metals and type of lattice. 
For example a-brass generally have below 35 atomic percent Zinc, a single phase and 
is made of a face-centred lattice. At compositions in a narrow range around 50 atomic 
percent Copper, 50 atomic percent Zinc, the atoms occupy the sites of a body-centred 



cubic lattice (as show in Figure 1.1) forming /3-brass. A body-centred cubic lattice 
consists of two interlocking simple cubic sublattices, where one site of one lattice is 
centred in the centre of a cube of the other lattice. This system undergoes a phase 
transition at the critical temperature, Tc, of about 740 K. The distribution of atoms 
on these sites is disordered above a temperature T^. Below there is ordering with 
atoms of each kind preferentially distributed on one of the two simple cubic sublattices 
of the body-centred cubic sublattices of the body-centred cubic lattice (/?' phase). 



1.4 Related works 

In this research we use a simulation program which uses a Metropolis-Hastings algo- 
rithm. This algorithm produces local changes in the system, which means the con- 
figuration is changed only in a small neighbourhood of each update. While there are 
advantages to this method, near the critical point there is a slowing down phenomenon. 
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Figure 1.1: Simplified representation of tlie body centred cubic of the /3' phase of 
brass where the red atom is Copper and the blue atoms are Zinc. 



This means that with an increase in lattice size the computation time increases rapidly 
so the numerical study has been restricted to small systems. To study larger systems 
it has been suggested by Bloete et al. in |BHL02j that cluster Monte Carlo methods 
can be used. In these methods, instead of flipping each individual site per iteration, 
clusters of sites are now flipped. This allows more efficient algorithms, such as those 
from Swendsen and Wang |SW87j . Bailhe and Coddington |BC91] and Wolff |Wol89j . 
The Wolff single-cluster algorithm |Wol89] stands out because of its simplicity with 
only one cluster being formed and flipped at a time. However, as cluster analysis is not 
as easy to generalise as the local interaction methods, its applicability range is limited. 
One possible pitfall of cluster algorithms is in the situation where the clusters tend 
to occupy practically the whole system, resulting in only trivial changes of the spin 
conflguration, and in a limited efficiency. 

The work of Wu and Lin has be extended by Strecka to mixed spin Ising models, in 
his papers |Str06j and |SvD06] . Mixed spin models are systems where the spins have 
different magnitudes. These mixed spin models have a much richer critical behaviour 
compared to the single spin models considered in this thesis. They are also useful as 



they describe the simplest ferrimagnets (see Section 2.5) and so have a wide range of 



practical applicability. In his solo paper |Str06j . Strecka investigates the spin- 1/2 and 
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spin-3/2 Ising model on the Union Jack lattice. In the article he presents a formulation 
of the system and its equivalency to the eight-vertex model, and then performs numer- 
ical simulations to investigate the critical points. In |SvD06j . he and his co-authors 
investigate the mixed spin-(l/2, S") Ising model, and compare the integer versus half- 
odd-integer spin-S* case. In the paper they extend Strecka's previous work and examine 
how critical exponents may depend on the quantum spin number. For example, the 
mixed spin-(l/2, 1) model exhibits quite different variations of the critical exponents 
when compared to the (1/2, 3/2) version. 



Chapter 2 

Theory Revision 



In this chapter we will present a brief review of some of the background theory that will 
be used in later sections. First we will look at the topic of thermodynamics, discussing 
both the motivation and the laws of thermodynamics. Then we will move on to look 
at some important topics in statistical physics like the partition function. Finally we 
will present an overview of the chemistry of magnetic substances. 

The topics in this section are commonly known results that form a basis for a great 
amount of work. As such discussions on these topics in statistical physics can be found 
in |Man88[lPB94] and for magnetochemistry |Ear68[ [HrcOSl ISha92| ISALOO] . along with 
many other books on both subjects. 

2.1 Motivation for thermodynamics 

A macroscopic system has many unmeasurable quantities that affect the overall be- 
haviour of the system. Thermodynamics concerns itself with the relation between a 
small number of variables that are sufficient to describe the bulk behaviour of the sys- 
tem in question. In the case of a magnetic solid the variables are the magnetic field 
B, the magnetisation M, and the temperature T. If the thermodynamic variables are 
independent of time, the system is said to be in steady state. If moreover there are 
no macroscopic currents in the system, such as a fiow of heat or particles through the 
material, the system is in equilibrium. 

A thermodynamical transformation or process is any change in the state variables 
of the system. A spontaneous process is one that takes place without any change in the 
external constraints on the system, due simply to the internal dynamics of the system. 
An adiabatic process is one in which no heat is exchanged between the system and its 
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surroundings. A process is isothermal if the temperature is held fixed. A reversible 
process follows a path in thermodynamic space that can be exactly reversed. If this is 
not possible, the process is irreversible. 

2.2 Laws of thermodynamics 

The laws of thermodynamics are relations that have been found to describe the prop- 
erties of a system when going through a thermodynamic process. These laws govern 

properties such as the direction that a process progresses in (which is stated by the 
second law). They also help to relate thermodynamics to mechanical ideas, such as the 
conservation of energy. 

Zeroth law 

Many quantities like pressure and volume have direct meaning in mechanical systems. 
However some basic quantities of statistical mechanics are quite foreign to mechanical 
systems. One of these concepts is that of temperature. Originally temperature relates 
to the sensations of hot and cold of a particular system. A remarkable feature of 
temperature is its tendency to equilibrium. One can see this when taking a soda can 
out of the fridge. When it is initially out of the fridge it is cooler than the surrounding 
room, but it slowly warms up to be in equilibrium, that is the same temperature, with 
the surroundings. This type of equilibrium is often referred to as thermal equilibrium. 

The zeroth law is used to assign a measurable value to this concept of temperature. 
Formally the law can be stated as: 

If a system A is in equilibrium with systems B and C, then B is in 
equilibrium with C. 

That is if two systems are in the same temperature as a third common system then 
they arc the same temperature as each other. In general thermometers are systems 
where a property that depends on its degree of hotness can be measured equally, 
such as the electric resistance of a platinum wire. Once calibrated with some known 
temperatures, such as the ice and steam points of water, and interpolating linearly 
for other temperatures, one can use this system to measure the temperature of other 
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systems. Of course temperature scales depend on the particular thermometer used. 
This arbitrariness can be removed by developing an absolute temperature scale using 



the second law of thermodynamics (2.2) 



First law 

The first law of thermodynamics restates the law of conservation of energy. However, 
it also partitions the change in energy of a system into two pieces, heat and work: 



dE = dQ-dW. (2.1) 



In (2.1) dE is the change in internal energy of the system, dQ the amount of 
heat added to the system, and dW the amount of work done by the system during an 
infinitesimal process. Apart from the partitioning of energy into two parts, the formula 
distinguishes between the infinitesimals dE and dQ, dW. The difference between the 
two measurable quantities dQ and dW is found to be the same for any process in which 
the system evolves between two given states, independently of the path. This indicates 
that dE is an exact differential or, equivalently, that the internal energy is a function 
based on the state of the system. The same is not true of the differentials dQ and dW, 
hence the difference in notation. 

Consider a system whose state can be specified by the values of a set of the sate 
variables Xj (for example the magnetisation, et cetera) and the temperature. Ther- 
modynamics exploits an analogy with mechanics and so, for the work done during an 
infinitesimal process, 

dW = -J2Xjdxj 

3 

where the Xj's can be though of as generalised forces and the X j S cLS generalised 
displacements. 



Second law 

The second law of thermodynamics introduces the entropy, S, as a state variable and 
states that for an infinitesimal reversible process at temperature T, the heat given to 
the system is 

dQrev = TdS (2.2) 
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while for an irreversible process 

dQ irrev 

< TdS 



If the only interest is in thermodynamic equilibrium states, (2.2) can be used and then 
the entropy 5* can be treated as the generalised displacement that is coupled to the 
force T. The above formulation of the second law is due to Gibbs. 

Two equivalent statements of the second law of thermodynamics are: The Kelvin 
version 

There exists no thermodynamic process whose sole effect is to extract a 
quantity of heat from a system and convert it entirely to work. 

and the equivalent statement of Clausius 

No process exists in which the sole effect is that heat flows from a 
reservoir at a given temperature to a reservoir at a higher temperature. 



Third law 

The third law of thermodynamics is of a more limited use, with its main applications 
being in chemistry and low temperature physics. It originated in Nernst's (1906) study 
of chemical reactions, and is sometimes called Nernst's Theorem. While it is of some 
interest in the report and it will briefly discussed, a more detailed treatment can be 
found in |Wil611 ItHM] . 

The third law of thermodynamics is concerned with the entropy of a system as the 
temperature approaches absolute zero. Stated simply it says that at absolute zero all 
bodies will have the same entropy. This means that at absolute zero a body will be in 
the only possible energy state, so would possess a definite energy (called the zero-point 
energy and would have zero entropy in this state. This value is left undetermined by 
purely thermodynamic relations such as the previous laws, which only give entropy 
differences. However a postulate related to the second law, although independent from 
it, is that it is impossible to cool a body to absolute zero by any finite process. Although 
you can get as close as you desire, you can not actually reach the absolute zero. 

An alternative version of the third law is given by Gilbert N. Lewis and Merle 
Randall in 1923 
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If the entropy of each element in some (perfect) crystalhne state be 
taken as zero at the absolute zero of temperature, every substance has a 
finite positive entropy; but at the absolute zero of temperature the entropy 
may become zero, and does so become in the case of perfect crystalline 
substances. 

Of course there are substances where at absolute zero there exists a residual entropy, 
such as carbon monoxide, but this is because there exist more than one ground state 
with the same zero point energy. 

2.3 Partition functions 

The partition function {Z) jMan88] is a quantity in Statistical Mechanics that encodes 
the statistical properties of a system in thermodynamic equilibrium. It is a function 
of temperature and other parameters, such as the volume enclosing a gas. Using the 
partition function of a system one can derive quantities such as the average energy, 
free energy, entropy or pressure from the function itself and its derivatives. 

There are many forms of the partition function depending on the type of statisti- 
cal ensemble used. The two major versions are the Canonical partition function and 
the Grand canonical partition function |PB94] . The canonical partition function ap- 
plies to a canonical ensemble, that is a system that is allowed to exchange heat with 
the environment at a fixed temperature, volume and number of particles. The grand 
canonical partition function applies to the grand canonical ensemble, where a system is 
allowed to exchange heat and particles with the environment, at a fixed temperature, 
volume and chemical potential. Of course for other circumstances one can define other 
partition functions. In this report we will be using the canonical partition function 
exclusively, so we will briefly develop this here, and refer to it from now on as the 
partition function. 

2.3.1 Definition 

Assume we are looking at a thermodynamically large system that is in constant thermal 
contact with its surroundings, which have temperature T, and the system's volume and 
number of particles have been fixed, in other words a system in the canonical ensemble. 
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If we label the microstates (exact states) that the system can take with s {s = 1,2,...) 
and the energy of these microstates as Eg then the partition function will be 

Z = ^e-^^^ (2.3) 

s 

where the "inverse temperature" /3 is defined as 

KbT 

In quantum mechanics, the partition function can be written more formally as the 
trace over the state space 

Z = Ti (e-^^) 
where H is the quantum Hamiltonian operator. 

2.3.2 Using the partition function 

It may not be obvious why the partition function is of such importance to Statistical 
Mechanics from the definition above. The ability to calculated the microstates of a 
system using a particular model, and from those to calculate a function that can be 
used to calculate the thermodynamical properties of the system is a very powerful 
tool. The partition function has a very important statistical meaning which allows this 
relation to the thermodynamical properties. The probability Ps that a system is in 
state s is 

Ps = ^e-^^^ (2.4) 

^-l3Es jg ^Yie well known Boltzmann factor or Boltzmann weight |Man88] . As we can 
see the partition function is playing the role of a normalising constant for this, which 
ensures that the probabilities add up to one: 

s s 

It is for this reason Z is called the "partition function"; it weights the different mi- 
crostates based on their individual energies. The reason for the partition function to 
be called Z is because it is the first letter of the German word for sum over states. 
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"Zustandssumme" . 

We can go on to use the partition function to calculate the average energy of the 
system. This is simply the expected value of the energy, which of course is the total of 
the energies of the microstates weighted by their relative probabilities. 

s s 

or equivalently, 



2.3.3 Entropy in statistical physics 

The idea of entropy was introduced above in the second law of thermodynamics, and 
it was even defined as a state variable, but what is entropy? The entropy of a system 
is a measure of the disorder of the system and is related to the energy of the system. 
Each microstate, which are the states at the microscopic level, have statistical weight 
related to the energy of the system. The entropy of these microstates can be defined 
from these statistical weights by using the following formula, 

S^kBlnn{E) (2.6) 

where Q{E) is the number of the microstates with energy E. The process in which the 
system is heated up without any work is reversible. Hence 

dE = dQrev = TdS. 



Therefore, by rearranging this formula 

dE T ^ ' 

where T is the temperature of the system. This shows that the entropy will increase 
slower as the temperature gets higher, but the entropy will never decrease as the 
temperature increases. 

As we have seen above, one can calculate many quantities with the partition func- 
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tion, and the entropy of the system is one of those quantities. 

S = -kB (In Z + /3 (E)) = ^ {ksT In Z) 

2.4 Phase transitions 

As we are studying the physics of phase transitions we should, at this point, define 
what we mean by a phase transition. To do this first let us define what we mean by a 
phase. 

A phase is a homogeneous part of a system bounded by surfaces across which 
the properties change discontinuously. In the most general situation each phase will 
contain several components, that is, it will contain several different species of molecules 
or ions. For example a gaseous phase might consist of a mixture of gases. To specify the 
properties of a phase we must then specify the concentrations of the various components 
in it. Then we have the possibility of transfer of matter between different phases. This 
matter has then gone through a phase transition. 

A first order phase transition is normally accompanied by the absorption or libera- 
tion of latent heat, while a second order phase transition has no latent heat. An example 
of a first order phase transition is between water and steam (gaseous water), and an 
example of a second order phase transition that exhibited by is the two = dimensional 
Ising model on a square lattice (magnetisation). There is a critical temperature and 
energy at which the properties of the two phases of the transition are identical. A 
phase transition is easily seen on a graph of the spontaneous magnetism versus the 
temperature of the system as a discontinuity of the data, or a point where the temper- 
ature remains constant while the energy input increases. Further information on phase 
transitions can be found in |Man88] . 

2.5 Magnetic substances 

Substances are made up of smaller particles that cannot be divided without loss of 
the description of the substance, that is, if these particles were divided, the parts 
would not exhibit all the properties of the substance. Depending on the bonding of 
the substance, these particles may contain smaller particles called ions or it may be a 
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molecule of the substance. Ions are of a similar configuration as atoms, with electrons 
orbiting around a nucleus, but have a different electron configuration. If the ion has 
more electrons than the base atom, then the ion will have a negative charge, and with 
less electrons will have a positive charge. Molecules are collections of atoms that share 
their electrons among each other by using overlapping orbits, which is called a covalent 
bond. As it is sufficient to describe the properties of the molecules or ions, because 
the substance is then in its most basic form, we shall concentrate on the description 
of these particles initially. Then we shall look at a crystal lattice structure of atoms 
or ions to investigate the possible properties that might be found as a result of the 
arrangement of the particles. It is sufficient to only describe a crystal lattice with ionic 
bonding as it can be generalised for any other structure with any type of bonding. 

For a molecule or ion to be paramagnetic |Ear68l IOrc03] . that is, attracted by a 
magnetic field, it must contain one or more unpaired electrons. Substances that are 
repelled by magnetic fields are called diamagnetic and contain no unpaired electrons. 
The magnetic moment of a paramagnetic species increases with the increase in the 
number of unpaired electrons. Liquid oxygen is easily shown to be paramagnetic by 
pouring it between the poles of a strong magnet, when it is attracted by the field and 
fills the gap between the poles. 

Substances in which the paramagnetic species are separated from one another by 
several diamagnetic species are said to be magnetically dilute |Sha92j . When the para- 
magnetic species are very close together (as in metals) or are separated only by an 
atom or monatomic ion that can transmit magnetic interactions they may interact 
with one another. This interaction may lead to ferromagnetism (in which large do- 
mains of magnetic dipoles are aligned in the same direction) or antiferromagnetism (in 
which neighbouring magnetic dipoles are aligned in opposite directions) |SAL90j . Fer- 
romagnetism leads to greatly enhanced paramagnetism as in iron at temperatures up 
to 1041 K (the Curie temperature), above which thermal energy overcomes the align- 
ment and normal paramagnetic behaviour occurs. Antiferromagnetism occurs below 
a certain temperature called the Neel temperature; as less thermal energy is available 
with decrease in temperature the paramagnetic susceptibility falls rapidly. More com- 
plex behaviour may occur if some moments are systematically aligned so as to oppose 
others, but to a finite magnetic moment: this is Ferrimagnetism. The relationship 
between paramagnetism, ferromagnetism, antiferromagnetism and ferrimagnetism is 
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illustrated in Figure |2.1 



(a) 



tititlt 



ititlt 



Figure 2.1: Pictorial representation of the relative magnetic spin directions of (a) 
paramagnetism, (b) ferromagnetism, (c) antiferromagnetism (d) ferrimagnetism. Here 
the relative magnitude of the spins is shown by the size of the arrow. Taken from 
Sharpe. |Sha92| 



2.6 Modelling magnetic substances 



Observing the properties of magnetic substances described in Section 2.5 is difficult as 



one has to examine the system on the quantum level with all the issues that brings. 
Instead, it is easier and cheaper to use a statistical model of the substance in a particular 
configuration^ such as the Ising model. In this section, we shall concentrate on the 
energy that the particular configurations have, as this is crucial to the spontaneous 
magnetisation of the system. 

The Ising model consists of a regular array of lattice sites in one, two or three 
dimensions, in which each site can be occupied in one of two ways. Each lattice site 
is occupied by an atom that can only exist in one of two spin states: +, with the spin 
in the direction of the magnetic field Bq, or — , against the field. The potential energy 
of a single dipole or spin is —itlBq if it is oriented with the field (+), and +mi?o if 
it is oriented against the field (— ), where m is the magnetic moment of an individual 
atom . Let the state of the j*^ lattice site be denoted by cxj, which equals +1 for a 
+ state and —1 for a — state. In terms of these cTj, the potential energy due to the 
external field is —niBQ^aj. To simplify this formula, for the later sections we shall 
use the following definition; h = uiBq. It is also assumed that there is an interaction 



"'^that is, a given temperature, lattice configuration and external magnetic field 
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eij between nearest neighbour atoms, which are in sites i and j, where i and j can 
be + or — . In terms of aj, can be written as —JaiOj where J is the interaction 
strength between spins. Note that: if the spins are parallel, aiaj is positive; if there 
are anti parallel, (TjCTj is negative. So the nature of the interaction is determined by 
the sign of J . If J > 0, parallel alignments are more stable and the model will describe 
ferromagnetism. If J < 0, an opposed alignment is the more stable and this will lead to 
antiferromagnetism. These interaction energies are similar to those in chemical bond 
theory |Sha92j . 

The total energy of a given configuration of N sites on a particular lattice type and 
variables, that is a given set of {cTj} is then 



where the first summation is over all nearest-neighbour pairs. This is the energy expres- 
sion that characterises the Ising model of a magnetic system. The canonical partition 
function is the summation over all configurations, weighted by exp {—E/ {kBT)) or 

ry ST- ^ I £'((Ji,Cr2, . . . ,CTAr) \ 

^ ■ ■ ^„ ' I 1^ i ' ' 

CTl=±l (Tjv=±l ^ ■' 

where fcs is the Boltzmann constant. The number of terms in this equation is 2^, 
because each of the sites can take on two values. 

In a certain sense, the Ising model is a simpler system than a non-ideal gas or 
a liquid, because the interacting particles are allowed to be situated only at discrete 
lattice sites. On the other hand, the model is difficult enough that the exact solution 
in three dimensions has yet to be found, although the two-dimensional model has been 
solved exactly in the absence of a magnetic field. 



N 



N 





Chapter 3 

The Ising Model 



In this chapter and the next we will look at the development of the solution of the Ising 
model for various types of lattices. We will first look at the one dimensional model, to 
show a simpler form of the model being exactly solved. Using this result we will then 
look at the two-dimensional model on a square lattice. Both the square lattice model 
and chain model are special in terms of Statistical Mechanics as they are examples of 
the few models that can be exactly solved. At the end of this chapter we will quickly 
cover the Ising model in higher dimensions and why these models may never be exactly 
solved. Also towards the end of the chapter we will quickly discuss other methods of 
solving and approximating solutions to the Ising model. 

3.1 Calculation for one-dimension 

To make the development easier we shall start with the one- dimensional model, or 
chain, with no external magnetic field. This will allow us to see how the interactions 
in the chain affect the properties. Later we will explain the model in an external field 
so as to improve the model for its use in applications. 

— e « o 

Figure 3.1: Pictorial representation of the one-dimensional Ising chain. Here the 
circles represent the atoms in the chain, and the lines represent the interparticle inter- 
actions with strength J. 

In our model, each configuration occurs with the probability proportional to e~^^, 
where E is the energy of the system. First, we shall work out the energy, or rather 
what we shall call the energy which is equal to the enthalpy as there will be no work 

21 
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done by the system, or each configuration which is given by the following Hamiltonian 
formula; 

= - J ^ o-jCTj+i (3.1) 

i=l 

where J is a constant. The term being summed is the relative spin configuration of 
the neighbouring elements. 

The canonical partition function [Zn) is given by 

j- iV-l >| ^ 

■■■ exp<l3jY \ where /3 = — — . (3.2) 



Zn — 

(Tl=ltl (T]V=il k j=l 



We can see from the equation that the spin of the last position appears only once in 
the equation and so we have the following formula irrespective of the value of crjv-i: 

Y^ exp {/^JcTAT-iCTAr} = 2 cosh/3J. (3.3) 

o-jv=±l 



Substituting this equation into (2.3) we get 

Z]si = [2 cosh /3 J] Zn_i 

Now iterating this operation we can reduce the equation to 

Zn = (2cosh/3J)^~^Z2 (3.4) 

where ^2 = exp {/3J(Ticr2} = 4 cosh /3J. (3.5) 

(Ti=itl a"2=itl 



So combining (3.4) and (3.5) we get the final formula for Z^: 

Zn = 2{2cosh(3jf~^ (3.6) 
The Gibbs free energy of the system is defined by 

G = -kB In Zn = -keT [In 2 + (A^ - 1) In (2 cosh (3 J)] 
In the thermodynamic limit (that is, the limit of the free energy as tends to infinity). 
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only the term proportional to is important and so the formula becomes 



G ^ -NkeT In {2 cosh 13 J). 



(3.7) 



These calculations give us the exact solution to the model while not in the presence 
of any magnetic field. We can adapt these formulae and find the solution for a model 
in an external magnetic field. However, one extra quantity we have to consider is the 
end effects of this chain when in this field. These end effects could make the equations 
more difficult than they would need to be, so we shall try to 'fix' the model so that we 
can ignore the end effects. We can safely 'fix' the end effects as they do not matter in 
the thermodynamic limit. This is because we assume the chain to be infinitely long in 
both directions without an end. 

One way to 'fix' the problem of the end effects is to consider the chain as a ring, 
that is, to set periodic boundary conditions. So we shall assume that the N*^ spin is 
connected to the first. Then the Hamiltonian becomes 



N 



N 



{3.1 



i=l 



where the spin labels run modulo (that is, N + i = i). This formula can then be 
simplified to 



N 



i=l 



h , , 



The partition function for this new Hamiltonian can be written as 



N 



or 



E ■■■ E /'E 

cri=±l iTjv=±l I, 1=1 



Jaiai+i + - ((7j + (7i+ij 



N 



<7i i=l 



(3.9) 



We can see that this formula has the potential of being very large, so at this point we 
introduce a 2 x 2 transfer matrix. 



Pu 
P-ii 



Pi-i 
P-i-i 



(3.10) 
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where the elements are defined as follows 



Pu 

P-i-i 
Pi-i = P-u 



oP{J-h) 



(3.11) 



We shall now use this transfer matrix to describe our partition function in terms of a 
product of these transfer matrices: 



7 = \ ^ P P ■ ■ ■ P = Tr P 



N 



(3.12) 



We can see that the matrix P can be diagonalised and the eigenvalues Ai and A2 are 
the roots of the determinant 

det (P - AI) = (3.13) 

where I is the 2x2 identity matrix. Similarly the matrix has the eigenvalues Xi 
and A^, and the trace of P^ is the sum of the eigenvalues: 



— Xi + A^. 



(3.14) 



The solution of (3.13) is 



Ai 2 = e^-^ cosh I3h ± \ e^'^'^ sinh^ l3h + e-'^P^ 



(3.15) 



We note that the eigenvalue associated with the positive root of equation (3.13), Ai 



is always larger in magnitude than the negative root eigenvalue. Now putting this 
partition function into the equation for the free energy we get 



G = -keT In (Af + A^) = -ksT \N\nXi + In 



1 + 



N 



(3.16) 



This approaches —iV/c^T In Ai as goes to infinity. Now taking this to the thermo- 
dynamic limit, 



G = -NknT In 



e^^ cosh (3h + Je^f^J sinh^ (5h + e-^^J 
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In the special case h = we obtain the resuh for the non-magnetic field case, so we 
have found the general formula. We may compute the magnetisation from 

1 dG kBTdXi 

'^ = ^'''^ = -Ndh=^m^- 

After some (straightforward) manipulations we find 

sinh/3/i 

m = — , (3.17) 
Vsinh^ (3h + e-^PJ 

We see that for h = there is no spontaneous magnetisation at any non-zero temper- 
ature. However, in the limit of low temperatures 



sinh^ I3h > e 



-4/3J 



for any h ^ and only a very small field is needed to produce near saturation of 
the magnetisation. The zero- field free energy will, in the limit as T approaches zero, 
approach the value G = —NJ corresponding to completely aligned spins. We could 
thus say that we have a phase transition at T = 0, while for T ^ the free energy is 
an analytic function of its variables. 

This is a very interesting result as phase transitions do not occur in the one- 
dimensional model at a finite temperature. If we look at the model with an initial 
state we should be able to intuitively illustrate why the phase transition is impossible. 
The model is set up so (Tj = 1 if i < / and cxj = — 1 if z > /, with the ground state 
Eq = —{N — I) J. There are N — I such states, all with the same energy E = Eq + 2 J. 
At temperature T the change in free energy will be 

AG = 2J -kBT\n{N -I). 

This quantity is less than zero for all T greater than zero in the limit where ap- 
proaches infinity. The expectation value of the magnetism of the system is zero by using 



our formula given in (3.17). As the system is at equilibrium if the value of magnetism 
is zero, the change in free energy would be zero. So we have a contradiction and there 
can not be a phase transition between the equilibrium phase and the ferromagnetic 
phase. 
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3.2 Square lattice two-dimensional calculation 



As we have mentioned in Section 1.2 of this chapter, the Norwegian American chemist 



Lars Onsager presented the first exact solution of the two-dimensional rectangular 
lattice Ising model without an external magnetic field |Uns44j . While this solution 
is historically important, it is also quite mathematically formidable. Since Onsager's 
solution there have been a number of more transparent solutions presented |Yan52t 
IKau49t IKW52j . and in this section we will look at the method of Schultz, Mattis and 
Lieb |SML64l IPB94j . This method has been chosen as it follows the same form as the 
calculation for the one-dimensional model shown above, and it will relate well to our 
later development of the two-dimensional model on the Union Jack lattice. 



■Q 
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Figure 3.2: Pictorial representation of the two-dimensional Ising model on the square 
lattice. Again the lines represent the interparticle interactions, with the horizontal 
interactions having strength Ji and the vertical interactions have strength J2. The 
sites are shown by circles. 



3.2.1 Transfer matrix 

In the calculation in the one-dimensional model we used a transfer matrix approach 
|NM53j . and this method can be easily extended to be used for the two-dimensional 
model. First let us look at the one-dimensional model and this time formulate the 



solution in a slightly different way. First let us rewrite equation (3.9) in the following 
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way 

where K = 13 J. 



Taking the two orthonormal states 1+1) and | — 1) we can form a basis and rewrite 
the Pauh operators as, 

^.-(' ' ] .^=f° '] .-f" ° I (3.19) 

with (Jx = cr+ + a~ and cxy = —i{a^ — a^). The Boltzmann weight exp{/3/icrj} can be 
expressed as a diagonal matrix, Vi, in this basis: 

(+1 |Vi| + 1) - (-1 |Vi| - 1) = e-^'^ 

or 

Vi = exp(^/i(7z). (3.20) 

Wc can also define the operator V2 corresponding to the nearest-neighbour coupling 
by its matrix elements in this basis: 

(+IIV2I + I) = (-l|V2|-l) = e^ 
(+IIV2I-I) = (-IIV2I + I) = e-^. 

Therefore, 

V2 = e^l + e-^ax = A{K) exp{K*ax} (3.21) 

where in the second step we have used the fact that {(Jx)^"' = 1- The constants A{K) 
and K* are determined from the equations 

A cosh K* = 

AsinhK* = e"-^ (3.22) 
or tanhX* = exp(— 2ir), A — a/2 sinh{2ir}. Prom these results we can rewrite the 
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partition function as: 

2" = ^ (/il |Vi| /i2) (/i2 IV2I /is) (/U3 |Vi| /i4) . . . (/i2Af IV2I /il) 

{^=+1,-1} 

= Tr (ViVa)^ = Tr (V2/'ViV2/')^ = Af + (3.23) 
where Ai and A2 are the two eigenvalues of the Hermitian operator 

V = iyy^YiWy^) = ^2 sinh(2i^)e^*"^/2e'^'*"^e^*"^/l (3 24) 

We have managed to get to this symmetric form of the transfer matrix V by using 
the property of the invariance of the trace of a product of matrices under a cychc 
permutation of factors. 

It is interesting to note that in this method we have taken a one-dimensional classi- 
cal statistics problem and transformed it into a zero-dimensional quantum-mechanical 
ground state problem. This is not the only time this transformation can be achieved 
and the result is quite general jSuz76l ISuz85j . There is a correspondence between 
the ground state of quantum Hamiltonians in — 1 dimensions and classical partition 
functions in d dimensions. This can sometimes be exploited, for example, in numerical 
simulations of quantum-statistical models. 

We can now generalise this method to the two-dimensional Ising model. This 
time consider an M x M square lattice with periodic boundary conditions and the 
Hamiltonian 

H = — J ar,cO'r+l,c ~ J Or^cf^r^c+l (3.25) 
r,c r,c 

where the label r refers to rows, c to columns, and a^+A/.c = c"r,c+M = CTr.c- Here the 
first sum contains only interactions in column c and the second sum is the coupling 
between neighbouring columns and will lead to a non-diagonal factor in the complete 
transfer matrix. 

Now as we did with the one- dimensional case, we introduce the 2^ basis states 

\^J) = 1/^2) • • • l/iAf) (3.26) 

with /ij = ±1 and M sets of Pauli operators {(Jjx-, o'jy, o'jz) which act on the j*^ state 
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in the product, that is, 

O'jZ 1/^1, /^2, • • • A'j, • • • /^m) = /ij 1/^1, /^2, • • • /^j • • • A^m) 
(T+ |/ii,/i2, . . ./ij, . . ./iM) = -1 l/il, Ai2, • • -/ij + 2 . . ./iA/) 

(^J \f^l,IJ'2, ■ ■ ■ fJ-j, ■ ■ ■ fJ'M) = S^.^i\fii,IJ,2, . . . fij -2 . . . flM) (3.27) 

Moreover, we impose the commutation relations [<Jja, o'mi3] = 0, a, (3 = X,Y, Z for 
j 7^ m. For j = m the usual Pauli matrix commutation relations apply. 

We can think of the index fij as the orientation of the jth spin in a given col- 
umn. From this we can see that the Boltzmann factors exp{K ar,c<^r+i,c} are given 
by the matrix elements of the operator Vi = exp{K o-jzo-(j+i)z}- Similarly, the 
matrix element 



l'^2| {yU'}) = ^/iM, /i Af-l, ••• ,Ail 

= exp{{M-2n)K} 



M 



n(e^l + e-^a,^) 



(3.28) 



where n of the indices {fi'} differ from the corresponding entries in {fi}. Then in a 
zero magnetic field the partition function of the two-dimensional model is given by 



{^J.l},{^^2},■..,{^^M} 



Tr (ViVs)^^ = Tr (vl^^ViVl^'^ 



M 



(3.29) 



Above the sum over each is over the entire set of 2 basis states. We can then 



write, using 3.21 and 3.22 



Va = (2sinh2ir)^/^exp 



(3.30) 



and as such the calculation of the partition function has been reduced to finding the 
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largest eigenvalue of the Hermitian operator 



V = vy^Yivy^ 



M \ ( M '\ ( M 



= (2 sinh 2K)^'^ exp | 5Z | | ^* ^^^^i+L^ J | ^* 5Z 

(3.31) 

This is still a non-trivial task since the factors do not commute with each other, and 
because the matrix V becomes infinite-dimensional in the thermodynamic limit. 



3.2.2 Transformation to an interacting fermion problem 

We now perform a rotation of the spin operators for future simplicity and let (jjz — > 
—UjXi <^jx <^jz for all j = 1, . . . , M. The eigenvalue are invariant under this rotation. 
Using (Tjz = 2o"^crJ — 1 and ajx = + '^7 '^^ arrive at 

K 



Vi = exp|irJ](a+ + a7)(a++^ + a7^i) 
V2 = (2sinh2ir)*'/'exp|2i^*^(^a;a7-il^ 



(3.32) 



Now following the method of Schultz et al |SML64] we perform a Jordan- Wigner 
transformation, which converts the Pauli operators into fermion operators (see |PB94j 
for a discussion of second quantisation). This step is necessary as some of the subse- 
quent canonical transformation are not possible for angular momentum operators. We 
write 

(T+ = exp <^ Tii c^.^c^ \ c] 

I m=l J 

aj = Cj exp <^ -TTi ^ aj^c^ I = exp <^ vri ^ cl^Cm > c (3.33) 

I m=l J I. m=l J 
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where the operators c, obey the commutation relations 



I '^i ' I _(_ 



I Cjf , Cm I _|_ 



0. 



The operator cJ^Cm, is the fermion number operator for site m with integer eigenvalues 
and 1. 

We can now express the operators Vi and V2 in terms of fermion operators using 



equation (3.33). The operator V2 presents no difficulties and is given by 



V2 = (2sinh2i^)*^/'exp i^2K* ^ (^c]cj - 



(3.34) 



In the case of Vi, there is a slight difficulty due to the periodic boundary condition. 
For j 7^ M we note that the term 



For the specific case j = M, 



{'^M + (^Ai) {(^t + 0^1 ) 



M-i "I r M-i 

exp <^ vri ^ c^cj > c\j (^c{ + Ci j + exp < Tri ^ c^j \ cm {c\ + ci 
j=i J I j=i 

M 



exp <^ ni c^jCj 



3=1 



exp vricjvjCM ( cjv/ + cm ) ( 4 + ci 



;-l)" (^cm - (^cj + ci 



where n = ^^cjcj is the total fermion number operator. The operator n commutes 

j 

with V2 but not with Vi. On the other hand, (—1)"" commutes with both Vi, and V2 
as the various terms in Vi change the total fermion number by or ±2. So we can 
write Vi in a universal way by considering the subspaces of even and odd total number 
of fermions, that is, 



{M 
K ^ [c] - (^cj+i + 9+1 



(3.35) 



32 



Chapter 3. The Ising Model 



(3.36) 



where 

Cm+1 = -ci, c[^+i = -c{ for n even 
cm+1 = ci, cj^^_^i = cj for n odd. 

This choice of boundary condition on the fermion creation and annihilation opera- 
tors allows us to recover translational invariance, so we can carry out the canonical 
transformation 



— Tc^- 



(3.37) 



with inverse 



-y 

Q 



(3.38) 



where we take q = jn/M to reproduce the boundary conditions in (3.36). Substituting 
this into the equations for V2 and Vi we get for n even, 



V2 = (2sinh2/s:)*^/^exp J 2/s:*^ (^aja, + at 



(2sinh2ir)*'/'J]V2, 

g>0 



and 



(3.39) 



Vi = exp < 2K l^cos q ^aj 



Qjq ~\~ Qj q(X — q 



i sin g ( a^alg + Oga.q 



g>0 



(3.40) 



In (3.39) and (3.40), the terms corresponding to q and —q have been combined. The 



resultant operators can be written as products and we can then see that bilinear op- 
erators with different wave vectors commute. This in turn allows great simplification 
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and the eigenvalues of the transfer matrix can be written as a product of eigenvalues 
of at most 4x4 matrices. For odd n we need the operators Vi^ and V25 for g = tt and 
q = 0. These are given by 

Vio = exp|2iraJ)ao} V20 = exp |2fs:* (^ajao - ^) } 
Vi^ = exp {-2Ka'la^} = exp {2K* (4a^ - |) } 

which are already in diagonal form and commute with each other. 



3.2.3 Calculation of eigenvalues 

We now go on to calculate the eigenvalues of the operator 

for q ^ and q ^ 71. As we are dealing with fermions, we only have four possible 
states: |0), |0), al^ |0) and ajal^ |0), where |0) is the zero particle state defined by 
Oq |0) = a_g |0) = 0. These states are eigenstates of V2. Since the operator Vi has 
non-zero off-diagonal matrix elements only when two states differ by two fermions, the 
problem reduces to finding the eigenvalues of in the basis |0) and |2) = a^aLq |0) . 
We note that 

Viguig |0) = exp {2K cos q} a^, |0) (3.42) 

and 

V^flO) = exp{-ir*}|0) 

V^f |2) = exp{X*}|2). (3.43) 
To obtain the matrix elements of Vi^ in the basis |0), |2), we let 



Vi,|0) = a(X)|0)+^(X)|2) 
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When we differentiate this with respect to K, we get 

da 



dK 

d/3 

dK 



= 2il3{K)smq 
= 4:P{K)cosq - 2ia{K)smq. 



(3.44) 



Subject to the boundary conditions a{0) — 1, /3(0) = 0, we solve these equations to 
find 



(0 |Vig| 0) = a(K)^ g2Kcosq (cosh2ir - sinh2ircosg) 
(2|Viq|0) = I3{K) = -ie'^^^°"^smh2Ksmq. 



(3.45) 



Using this method the matrix elements for (2 |Viq| 2) and (0 |Viq| 2) = (2 \Viq\ 0)* and 
obtain the matrix 



cosh 2K — sinh 2K cos q i sinh 2i^ sin q 

—i sinh 2i^ sin q cosh 2i^ + sinh 2K cos g 



(3.46) 



and 



exp{-K*} 
ex.p{K*} 



[V 



19J 



exp{-ir*} 
ex.p{K*} 



(3.47) 



The eigenvalues of this matrix are easily determined and can be given in the form 



= exp {2i^ cos q± e{q)} 



where after a bit of algebra e{q) can be defined by 



(3.48) 



cosh e(g) = cosh 2K cosh 2K* + cos q sinh 2K sinh 2K* . 

By convention we choose e(g) > 0. It can be seen that the minimum of the right hand 
side of this equation occurs as 5 — > tt and that for all q 



e{q) > e^nin = lim e(g) ^2\K-K* 



(3.49) 
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and also 



\ime{q) = 2{K + K*). (3.50) 



Now wc combine this information and consider the subspace on an even number of 
fermions. The wave vectors in this case do not include g = or g = tt, and it can be 
seen that the largest eigenvalue of for each g is A+. Thus the largest eigenvalue in 
this subspace is given by 



Ae = (2sinh2X)'^/' JJa+ 



q>Q 



(2 sinh 2K)^/^ exp i ^ [2 cos g + e(q)] 



= (2sinh2ii')'^/'exp|^ J]e(g)| (3.51) 

where we have used that X^gCOsg = and have extended the summation over the 
entire range, — tt < q < tt. 

It is slightly more difficult to work with the odd subspace. As before for g 7^ and 
g 7^ TT the maximum possible eigenvalue is A^. The corresponding eigenstates are all 
states with (—1)" = —1. So to make the overall state have the property (—1)" = — 1, 
we occupy the q = state and leave the q = n state empty, and as such obtain a 
contribution of (2 sinh 2K)^^'^ exp {2K} to the eigenvalue Ag. So the largest eigenvalue 
in the odd subspace is 

Ao = (2 sinh 2K)^/^ exp J 2X + ^ ^ e(g) I . (3.52) 

With further consideration of these results, it can be shown that Ag and Ag are de- 
generate. The critical temperature of the two-dimensional Ising model is given by the 
equation K = K* or by expression 

sinh = 1 (3.53) 

Kb^c 

which is approximately ksTc/ J = 2.269185 . . . 

The degeneracy of the two largest eigenvalues is neghgible in the dimensionless free 
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energy and only contributes an additive term of In 2. So for any temperature the free 
energy is given by 



^-^^ = /3,(0,T) = -iln(2sinh2ir)-^5^e(g) 

g 

- In (2 sinh 2K) - — f dq e(g) (3.54) 
2 47r 



2 

where the sum over the wave vectors has been converted to an integral. 



3.2.4 Spontaneous magnetisation 

While other thermodynamic functions can be derived, now we have the free energy of 
the system, which can be simplified with a bit more algebra. A full derivation of this 
result along with the derivation of the specific heat capacity of the model can be found 
in |PB94] . It is also possible to extend the theory discussed in this section to calculate 
the spontaneous magnetisation, as in |SML64j . The result is 



mo(r) 



d 

lim -7rrg(h, T) 



'1 



-tanh^^j)'^ 
letanh^^J 





(3.55) 



As T — )■ Tc, the limiting form of the spontaneous magnetisation is given by 



mo(T)^(T,-T)^/« = (Te-T)^ 



(3.56) 



At the critical point, the order parameter has a power law singularity. 



3.3 Triangular lattice two-dimensional model 

Now we add diagonals in one direction to the square lattice to obtain the triangular 



lattice seen in Figure [3^ This lattice configuration will be used as our "test" lattice in 
the later chapters. In this section we will state the results for the two site interaction 
spontaneous magnetisation and then go on to look at the three site interaction result. 
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Figure 3.3: Pictorial representation of the two-dimensional Ising model on the Trian- 
gular lattice. Interparticle interactions are shown here with lines, and have interactions 
strength Ji horizontally, J2 vertically and J diagonally. Sites are shown with circles. 



3.3.1 Two-site interactions on a triangular lattice 



Stephenson in his paper }Ste64j shows the development of the two site interaction 
result. As this development is similar to the development for the square lattice result, 



here we will state the final result for the spontaneous magnetisation in (3.57): 



M 




T <T, 

T>Tr 



(3.57) 



where 



[(1 



16 (1 + V1V2V3) {Vi + V2V3) {V2 + V3V1) (f3 + V1V2) ' 

Vi = tanh pji, V2 = tanh PJ2, V3 = tanh /3J. 



(3.58) 



Similar derivations can be seen in |Pot52l IGre62] . 



3.3.2 Three-site triplet interactions on a triangular lattice 



Having a result for the two-site nearest neighbour interactions, we now move on to 
look at the three site triplet interactions. When we consider these interactions the 
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Hamiltonian (3.25) becomes 





(3.59) 



where the first summation is over all 2N three-site triplet interactions, the second is 
over all 3A^ two-site nearest neighbour interactions and the third is over all N sites of 
the lattice. When any of the two of B, J, J' are zero, the free energy of the system 
can be evaluated exactly. In the case where B = J = we obtain the pure three-spin 
model, whose free energy was obtained by Baxter and Wu in |B W73j . 

In his approach to calculate the spontaneous magnetisation of the three-site inter- 
actions, Baxter |Bax75j considers the case when H = J' = 0. This allows a return to 
the normal two-site interaction triangular Ising model. For this he introduces a new 
parameter for the three-site spontaneous magnetisation, M3, and defines T^oo to be 
the ratio between the three-site correlator and the magnetisation in the infinite lattice 
case, i.e. 

7? = - 

Here the three-site correlator is the interaction between a triplet of sites ax, (Jy Oz on 
a triangular face, {a^ cXy cr^). From (62) of |Bax75j we have that 



This then means we can evaluate the spontaneous magnetisation of the three-site in- 
teractions as 



3.3.3 Three-site triplet interactions on a square lattice 

Now consider the case when J = V3 = 0. By removing the diagonal interactions we 
have the square lattice. M3 now becomes the three-spin magnetisation around a corner 
of the square lattice. Using the equations above we can obtain 



■00 



- + t;/ + f 2 + f 2 ^ + f 3 + f3 ^) 

-7: {viV2V3)~'^ [(1 + V1V2V3) {Vi + V2V3) {V2 + V3V1) {V3 + ViV2)]^ ■ (3.60) 




(3.61) 



Ms = M 1 - 



4g-4/3Ji-4/3J2 



(3.62) 
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where M is the spontaneous magnetisation of the square lattice. Other three-spin 
magnetisations have been evaluated by Pink jPin68j . 



Chapter 4 

The Union Jack Lattice 



Having discussed the two-dimensional model both on the square and triangular lattices 
in the last chapter, in this chapter we will focus on the Union Jack lattice or centred 
squared lattice. We obtain the Union Jack lattice by adding alternate diagonals to the 



squares of the square lattice, as shown in Figure 4.1[ Here the the nearest neighbour 




Figure 4.1: Pictorial representation of the two dimensional Ising model on the Union 
Jack lattice. Here, as before, the interparticle interactions are shown with lines. The 
red circles represent r-sites which interact with four particles. The f)hic circles represent 
cj-sites with eight interparicle interactions. The various interaction strengths are shown 
on the different interactions. 



interaction strengths can be one of six values Ji, J2, J3, J4, J, J' . The Hamiltonian 
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for this lattice is 



H 




(4.1) 



It is clear to see that this lattice (Ci) can be made of two sublattices, indicated by red 
and blue circles. The sublattice C" contains the r sites, red circles, with four nearest 
neighbour interactions and the sublattice C contains the a sites, blue circles, with 
eight nearest neighbour interactions. The partition function on £i is 



Vaks, Larkin and Ovchinnikov |VL066j first considered the Union Jack lattice Ising 
model as a system exhibiting a re-entrant transition in the presence of competing 
interactions. They considered a system where the interactions were symmetric, and 
obtained its free energy and a sublattice two-spin correlation function. This was later 
generalised by Sacco and Wu |SW75j as a 32-vertex model on a triangular lattice. 

Wu and Lin presented a simple formulation for the general Union Jack lattice model 
|WL87l IWL89] and showed that it is equivalent to a free fermion model |FW70] . Using 
this equivalence they obtained the free energy and the eight site interaction sublattice 
spontaneous magnetisation. This was extended for symmetric interactions by Lin and 
Wang |LW87j to produce a result for four interaction sublattice. Wu and Lin in their 
1989 paper |WL89] then extend the result further to produce results for both sublattices 
in the general interaction model. 



Although we will be presenting the general Union Jack lattice result from Wu and Lin 
|WL87| IWL89j in this chapter, it is useful to look at the results from the paper by 
Valks et al |VL066j . In their paper, they use a symmetric interaction model where 
the Jri r = 1, 2, 3, 4 are equal and J = J'. They showed that the phase of the system 




(4.2) 



4.1 The Valks et al result 



4.1 The Valks et al result 
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on the o"-sublattice can be determined by two parameters ai and 0:2, given by (4.3). 



ai 



(1 + 

e~^^ li - e 



a2 



2K _ g-2i^) 
-2K\ 



(cosh AKi + e- 



(4.3) 



where K = J/ksT and Ki = Ji/ksT. The critical temperature Tc is determined 
when «! = 1, which has solutions when —\Ji\<J. A second critical temperature 
T* is determined by 02 = —1, which has solutions when J2 < — 0.907 |Ji|. This is 



shown in Figure 4.2 With further analysis of the formulas given by Valks et al one 



1 I I I -p-i — ! — rj—i — 

\ ^ 

' \t^ / 

\C / 


1 1 1 — 1 — 1 — 1 — r~T — 1 — 
-1. ^.^-^ - 




-0.5 




«■ 


1 1 : 1 1 / i 1 > • 


— 1 — 1 — 1 — 1 — 1 — 1 ... 1 1 1 ... 



-1.0 



-0.5 



0.0 0.5 IjO „ 

q2J2/ccq,Jp2+cq2J2D2i'/2 



Figure 4.2: Phase diagram of the Vaks tt al. result on the Union Jack lattice Ising 
model showing the graphs of the critical temperature Tc, disorder point To, and second 
critical temperature T* . Below Tc the system is in a ferromagnetic phase. Below T* 
the system is in an ordered antiferromagnetic phase. Above these critical temperatures 
the system is in a disorder phase. Taken from |Ste70] . 

can classify the phases of the system. When — | Ji| < J2 < 0, there is an ordered 
ferromagnetic phase below Tc. Above Tc, there is disordered ferromagnetic phase up to 
a disorder temperature To, determined by ai = —02- At temperatures above T^, there 
is a disordered antiferromagnetic phase. This region extends to T = 00, if J satisfies 
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— 0.907|Ji| < J < 0. However, if — |Ji| < J < — 0.907 |Ji|, there is an intervening 
ordered antiferromagnetic phase between the lower and upper critical temperatures 

c • 



4.2 Computation on the Union Jack lattice 



Now we move on to look at the exact result for the general anisotropic Ising model 
on the Union Jack lattice from Wu and Lin |WL87t IWL89j . As the development is 
similar to that of the two-dimensional square lattice model shown in section 3^, we 
will present the result only. 



By summing over all the N/2 r sites on C" we can rewrite the partition function 
Zi as limit variables 



Ti=ibl Tjv=il V.i,j,k,l 



(4.4) 



where the product is over all N/2 faces of the sublattice each surrounded by four 
spins i, j, k, I. Here the Boltzmann factor u}{ai, o"2, cts, o'^) is obtained by dividing each 
of the diagonal interactions — J and — J' into halves, one for each of the two adjacent 
faces. This leads to 



u (a, b, c,d) = 2 exp 



(3J{ah + cd) (3J'{ad + hc) 



2 2 

This equation with its four spin variables gives us sixteen separate states which by 



cosh (a/3 Ji + h(3J2 + c/3 J3 + J4) . 

(4.5) 
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symmetry can be reduced to eight distinct expressions: 





= (jj 


(+ + ++) 


= 2e^-^+^'^' cosh (/3 ( Ji + J2 + J3 + ^4)) 


UJ2 


= OJ 


(+-+-) 


= 2e-^^-^-^' cosh (/3 ( Ji - J2 + - ^4)) 


UJ3 


= OJ 


(+--+) 


= 2e-^-^+^^' cosh (/3 ( Ji - J2 - ^3 + ^4)) 




= OJ 


(+ + — ) 


= 2e'^-^-^-'' cosh (/3 ( Ji + J2 - ^3 - ^4)) 


W5 


= OJ 


(+ - ++) 


= 2cOsh(/3(Ji- J2 + ^3 + ^4)) 




= OJ 


(+ + +-) 


= 2cosh(/3(Ji + J2 + J3- J4)) 




= OJ 


(+ + -+) 


= 2cosh(/3(Ji + J2-^3 + ^4)) 


UJs 


= OJ 


(- + ++) 


= 2 cosh (/3 (- Ji + J2 + J3 + ^4)) ■ 



(4.6) 



There can only be eight possible arrangements of the four spin variables and so we can 



consider it to be an eight-vertex model with weights (4.5). This eight-vertex model 



was proved to satisfy the free fermion condition by Fan and Wu |FW70] . 

OJ1OJ2 + OJ'iOJ4^ = WsWg + ^7^8 

and as such is a free fermion model with the partition function: 



(4.7) 



Z = -Zi. 
2 

Here the factor 1/2 takes account of the two-to-one mapping of the spin and vertex 
configurations. The spontaneous magnetisation of a free fermion model was given by 
Baxter in |Bax86j and is 



0, n-^ < 1, 



1 - 



7i 72 73 74 

I6U5 Uq OJj CJg 



(4.8) 
(4.9) 
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where 

71 = -tUi + U2 + UJ3 + 

72 = COi - UJ2 + UJ3 + 1^4 

73 = OOi + 002 - 003 + 

74 = OJl + U2 + UJ3 - 1^4: (4-10) 

The system exhibits an Ising transition at the critical point (s) 

= 1 (4.11) 

or equivalently, 

ui + U2 + UJ3 + U4 = 2 max {uji,uj2, W3, ^4} . (4-12) 



We can now go on to finding the spontaneous magnetisation for the other sublattice, 
C. While Lin and Wang presented a result for the symmetric case in 1988 |LW88j . here 
we will quickly review the result for the general case given by Wu and Lin in |WL89j . 
Their approach uses the three site triplet interactions and again, as in section |3.3.2 
we will simply present the result here. The r-sublattice magnetisation is given by 



(r) = (a) [A,234{K){F+ + F_) + A23,i{K){F+ - 



(4.13) 



We can see that this equation has a similar form to (3.61), and can be seen to be a 



multiple of the cr-sublattice value. In (|4.13|) 

Ai23,{K) -- 
^34l(i^) = 



sinh2(/3Ji + /3J3) 



^2G. iPJ) sinh 2/3 Ji sinh 2(3 J3 

smh2{(3J2 + (3Ji) 
^2G_ {/3J) sinh 2/3 J2 sinh 2/3 



and 



G_(/3J) = cosh2(/3Ji + /3J3) + cosh2(/3J2 - /3J4). 
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The calculation for and F_ is a little more involved. We start by calculating 



F± = W — (4.14) 



where 



A = 2UJ5UJQUJ7UJS {col +^2+ ^l^T) ~ {^1^2 + (^3i^4) {(^1(^3 + (^2^04) {C0l(^4 + (^2^03) 

B — U^UJqUJyUJs {oJ^UqUiUs — OJ1U2UJ3U4) 

C = {ul + ul + ujI + Ujf) ^ - 4 {u^Uq - ^7^^)^ 

D = ((^1 + cjI) {2u^uqU'jU^ — uJiU2iJJ2,U4) — UJ5UJQUJ7UJS {coj + uf) 

E ^ oof- 

We can relate and F_ with the following formula, allowing us to get values for each 
variable, 

F^F_ = (4.15) 

CJ\ljJ2 

We can compute the overall nearest neighbour magnetisation by taking the mean 
of the two sublattice magnetisations 

Mo=^((a) + (r)). 
4.2.1 Classification of phases 

At low temperatures the phase of the cr-sublattice can be classified based on the fol- 
lowing energy value. 



Ex = J + J' + I Ji + J2 + J3 + J. 

-E2 = -J- J' + |Ji- J2 + J3- J. 



-E4 = J- J' + |Ji + J2- J3- J4I (4.16) 
The sublattice is in a ferromagnetic phase when 



E\ < E2, E^, E4; 
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is antiferromagnetic when 

E2 < El, E^, E4; 

and finally is metamagnetic when 

E^ < El, E2, E4 or E4 < El, E2, E^. 

As the temperature rises, depending on the relative strengths of the interactions Jr, 
the occurrence of phase change is signified by one or more of the following equations 
being realised 

Wi + W2 + + 1^4 = 2 max {ui, U2, W3, ^4} . (4.17) 
A re-entrant transition occurs if any one equation admits two solutions. 



4.3 Theoretical analysis 

In this section we will plot the theoretical predictions of Wu and Lin, and compare 
those against the critical temperatures from the results of Stephenson. To do this we 
will plot results from each of the possible systems, along with a plot of the 7 from 



(4.9). Wu and Lin's predictions will be plotted against our simulation results later 
in Chapter |8j Here we are only concerned with identifying systems that need further 
investigation, due to conflict between the predictions of Vaks et al. and, Wu and Lin. 

For our first system we will look at the isotropic ferromagnetic case. This is the 
simplest system as all the interactions are the same, J„ = lOO/c^. The graph of this 



system is shown in Figure 4^: We can see that the prediction of this system moves 
from having a non-zero average magnetic spin to a zero average magnetic spin. The 
system has this phase transition at the critical temperature, around 400 Kelvin. 



This prediction agrees with the value obtained from ai in |4.3[ and a solution to 02 
does not exist. Intuitively this is the sort of graph we would expect. Although for this 
system both theories agree, we will go on to examine the changes in the gamma terms. 



This will set the groundwork for our later analysis. This plot is shown in Figure 4.4 
Here we can see that at low temperatures 71 is the only negative term, and that it 
equals zero at and then moves to being positive and converges to the level of the 
other 7. 



4.3 Theoretical analysis 
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Figure 4.3: Theoretical predictions of Wu and Lin for a ferromagnetic system with 
interaction strengths Jn = J = J' = WOks- Here the prediction for the cr-sublattice 



(eqn. (4.8)) is shown in black, the prediction for the r-sublattice (eqn. (4.13)) is shown 
in red and the prediction for the overall lattice is shown in blue. A phase transition 
can be seen just below 400 Kelvin. 




Figure 4.4: Graph of the separate 7 terms from (4.10) and the cj-sublattice prediction 
for the ferromagnetic system with interactions strengths Jn = J = J' = lOOks- Here 
73 is overlaid on 74 as they are equal. Note that at low temperatures 71 is negative 
and the other terms are positive. At the critical temperature 71 = 0. 
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Now we will move on to looking at a system that starts in a metamagnetic phase. 
For this system we need to set the diagonals to be opposite polarities, that is J = — J' = 
looks- Due to this the Vaks et al. result does not apply to this system. We also set 
the horizontal and vertical interactions to be J„ = 10k b- The plots for this system are 



shown in Figure 4.5 As a metamagnetic system only has an average magnetic spin when 



PrGdiction sigma 
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Figure 4.5: Results for a metamagnetic system with interaction strengths J„ = IO/cb, 
J = — J' = IOO/cb. (a) shows the theoretical predictions of Wu and Lin. We see that the 
result for the r-sublattice is zero across the temperature range while the a-sublattice 
shows a transition at just above 200 Kelvin, (b) shows the values of the separate 7 
terms. Here we see that 74 is negative at low temperatures and crosses the x-axis at 
the transition temperature. In this graph 71 is overlaid on 72 as they are equal. 



in the presence of an external magnetic field, we can see that something is incorrect 



here. From Figure 4.5b we can see that 74 is negative at low temperatures and behaves 
similarly to 71 in the ferromagnetic case. As 



74 = wi + W2 + ^3 - a;4 < 0, 



we note that uj^^ = uj{+ H ) is the dominant term. As this term would suggest an 

antiferromagnetic phase we would intuitively expect (cr) = 0. 

The system we next consider is that on an anisotropic antiferomagnetic system. 
For this system we set the vertical and horizontal interactions to be J„ = lOOfc^ and 
the diagonal interactions to be J = J' = — IOO/eb. As the system is antiferromagnetic 
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we expect it to have a zero average spin across the temperature range. The plots 



are shown in Figure 4.6 below. From Figure 4.6a we can see that the (r) prediction 
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Figure 4.6: Results for the antiferromagnetic system with interactions Jn = lOOfc^, 
J = J' = — lOOfce- (a) shows the predictions from the results of Wu and Lin. Note 
that the graph for the r-sublattice is above +1 which is physically impossible in our 
model. There is a transition below 200 Kelvin in all predictions, (b) shows the separate 
7 values for the system. Here 72 is identified as the dominant term for the system, and 
is negative while the predictions are non-zero. Again in this plot, 73 is overlaid on 74. 



produces physically impossible results. In our system the values of the spin can only be 
±1 and so a result greater than one is impossible. As the value of (r) is determined as 
a multiple of the value of (a) we should also look at this plot. We note that although 
our system is antiferromagnetic the prediction graph shows a non zero value for the 



magnetisation. As before, when we plot the individual 7 in Figure 4.6b we see that 
this time it is 72 that is negative before the critical temperature. This implies that 

(jj2 = uj{-\ 1 — ) dominates, and so again should have average spin of zero. When 

we compare these predictions to those of Vaks et al. we can see that the critical 
temperature here is T^, that is the disorder temperature where the system goes from 
being an ordered antiferromagnetic phase to a disordered antiferromagnetic phase. 

As one of the motivations for studying the Union Jack lattice is the property of 
re-entrant phase transitions, we will now move on to looking at one such system. The 
system we will now study is the anisotropic ferromagnetic case. For this type of model 
we will look at the system where the square lattice has interactions with strength 
Jn = IOO/cb while the diagonals will have interaction strengths of J = J' = —92^;^. 
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When we calculate the possible critical temperatures from Vaks et al. , we see that 
all three are possible, and there exists a re-entrant phase transition. The graph of Wu 



and Lin's prediction is shown in Figure 4.7 
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Figure 4.7: Results for the anisotropic ferromagnetic system with interaction 
strengths J„ = lOOA:^, J = J' = —92kB- This means that the horizontal and vertical 
interactions are equal, (a) shows the predictions from Wu and Lin for this system. 
Note after the first phase transition there is a second area of non-zero magnetisation. 
Here the prediction for the r-sublattice again has physically impossible results given 
the conditions of our model, (b) shows graphs of the separate 7 terms. Here we see 
that initially the 71 term is negative up to the first critical temperature. In the range 
of the second area of non-zero magnetisation we see that 72 is negative. As 73 is equal 
to 74 their plots are overlaid in this graph. 



Looking at Figure [4.7a| we see that in the set of curves between 45 and 100 Kelvin 
we have physically impossible results again for the (r) variable. Up to the first critical 
temperature we can see that the system behaves as we would expect a ferromagnetic 
system to, and then between the second and third critical temperatures we have the re- 



entrant phase. Examining Figure |4. 7b | we see that the first of these curves is determined 
when 7i is negative, suggesting a ferromagnetic system. During the re-entrant phase 
however, it is 72 which is negative, suggesting that its is an ordered antiferromagnetic 
region. In the paper of Vaks et al. |VL066j the re-entrant phase transition is defined 
to be between a disordered antiferromagnetic phase and an ordered antiferromagnetic 
phase. From only looking at the average magnetisation of the system it is not possible 
to see the transition from a disordered phase to an ordered phase. However the critical 



temperatures here are consistent with those predicted from (4.3). 
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However, the results from some anisotropic ferromagnetic systems do have agree- 
ment between both theories. To show an example of this, we will look at the system 
where the interactions on the square lattice are Jn = — lOO/cs and the diagonal inter- 



actions are J = J' = 92kB. This is illustrated below in Figure 4.8 
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Figure 4.8: Results for the anisotropic ferromagnetic system with interactions 
strengths Jn = —lOOks, J = J' = 92 A;^. These are the same interaction strengths 
as the last example but with the opposite sign, (a) shows the predictions of Wu and 
Lin for this system. Note that the overall system is antiferromagnetic with a zero aver- 
age magnetisation but the sublattices have spins of different signs. As the temperature 
rises, due to the shape of the different curves, an overall magnetisation increases up to 
the critical temperature of the phase transition, (b) shows the graphs of the individual 
7 terms. Here we see that only 71 is negative in the temperature range, becoming 
positive at the critical temperature. In this graph 73 and 74 are equal and as such 
their plots are overlaid. 



Looking at Figure 4.8a we can see that the curves are of opposite signs, but both 
display a ferromagnetic shape. Using the calculations from Vaks et al. we would expect 
the critical temperature to be at the value it is, and of the type it is. From examining 



Figure 4.8b we can see that it is only 71 which is negative at low temperatures. Equally 

as uji = u{+ + ++) = ), the result for (r) is acceptable, even though the 

overall system is antiferromagnetic. 

However there are some anisotropic systems, with non-uniform square lattice inter- 
actions where we see some interesting results. We shall now go on to look at a system 
where the horizontal interactions are Ji = J3 = lOOfc^/O.Q^, the vertical interactions 
are J2 = Ja = 100A;b/0.9 and the diagonal interactions are J = J' = IOO/cb. The graph 
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of this prediction is shown in Figure |4.9| below. 
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Figure 4.9: Results of the anisotropic ferromagnetic system with interactions hori- 
zontally of Ji = J3 = 100A;b/0.9^, vertically of J2 = J4 = lOOfcs/0.9 and diagonally 
of J = J' = woks- (a) shows the predictions of Wu and Lin for the system. Note 
that the predictions are all physically plausible being less than or equal to +1. In 
the r-sublattice prediction it can be seen that there is a little dip in the value before 
returning to the curve similar to the prediction function of the cr-sublattice. (b) shows 
the separate 7 terms. As with previous ferromagnetic systems 71 is the dominant term, 
and the graphs of 73 and 74 are overlaid. 



As we can see from Figure [4l9a the curves for r and subsequently the mean predic- 
tions deviate from the shape of the a curve at temperatures below 200 Kelvin. Above 
this temperature, the graphs show the expected ferromagnetic curves. As there are 
no values over +1 this is physically plausible. When we look at the separated 7 in 



Figure 4.9b the only term that crosses the x-axis is 71. This again suggests a purely 
ferromagnetic system, and shows the critical temperature where we would expect. We 
now look at a similar system where we swap around the interaction strengths of the 



vertical and horizontal interactions. The graph of this system is shown in Figure 4.10 



From Figure |4.10a we can see that again, the r and so the mean magnetic spin 
results are physically impossible given our system. However we can see from this figure 
and Figure |4.10b that the prediction for a remains the same. Intuitively we would 



expect this, as we have only rotated the system from Figure 4.9 by 90 degrees. Again 



we see that the predictions become more like we would expect above 200 Kelvin. 



4.4 Summary 



55 



< 0.6 



0-2; 



Prediction sigma 
Prediction tau 
Prediction mean 



°J 100 200 300 400 500 600 70D 800 
Temperature 




gamma 1 
gamma 2 
gamma 3 

lamma 4 

recJiction sigma 



(a) 



(b) 



Figure 4.10: Results for the anisotropic ferromagnetic system with interaction 
strengths Ji = J3 = IOO/cb/O.D, J2 = J4 = lOO^ij/O.g^ and J = J' = WOks- Note 



these interactions are those of Figure 4.9 rotated by 90 degrees, (a) shows the predic- 



tions of Wu and Lin for the system. Here we see that the prediction for the r-sublattice 
is physically impossible in our model below around 250 Kelvin, (b) shows the graphs 



of the separate 7 terms. Note that the result is identical to Figure 4.9b 



4.4 Summary 

We have seen in this section that there are a number of systems for which the predicted 
resuhs from Vaks et al. and Wu and Lin disagree. Equally we note that the Wu and Lin 
prediction only works when either 71 < or 71 72 73 74 > 0. These added conditions 
were not stated in either Wu and Lin's 1987 |WL87j or 1989 jWL89j papers. This 
suggests that these results need to be investigated further with other methods, such as 
mean field theory in Chapter [5] and numerical simulations in Chapter [8] 



Chapter 5 

Mean Field Theory 



In this chapter we will briefly discuss alternative ways of investigating the Ising model 
through the use of approximations. While exact solutions are useful for investigating 
systems under special conditions, such as in the absence of a magnetic field, approx- 
imations allow us to look at systems where an exact solution may not be possible. 
After this discussion we will look at one particular method of approximation, mean 
field theory, in more detail. Initially we will look at how mean field theory is defined 
for the triangular lattice model |Bax89j . which is a specific case of the mean field the- 
ory method for the isotropic lattice. After this we will develop our own equations for 
the Union Jack lattice. At the end of the chapter we will the perform some numerical 
analysis to investigate how well the approximation models our systems. 



5.1 Approximations 

Mean field theory is an older approach to the subject of phase transitions and generally 
give us a qualitative description of the phenomena of interest. There are many different 
approaches that can be use with mean field theory, but common to all is the identifi- 
cation of an order parameter. For this section we will be using the Curie- Weiss model 
|PB94] . In this method we approximate a system where nearest neighbours interact 
with each other, to a system where the particles do not interact with each other but 
instead with the mean spin of the entire lattice. This means that the approximation 
is affected less by the dimensionality or lattice configuration of the system and so we 



are able to avoid the issues discussed in Section 1.2 (work by jBar82t llstOOj ). This 
will allow us to investigate systems with external magnetic field in two dimensions. 
Alternatively, a similar approach where the free energy is expressed, then minimised 
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with respect to the order parameter, can be used. An example of this free energy 
minimisation approach is the Bragg- Wilhams approximation |BW34t lYewTOj . 

These methods can be improved to gain numerical values for the properties of the 
Ising model, by such approximations as those presented by Bethe in |Bet35j and its 
extension by Guggenheim |FG40] . However, the asymptotic critical behaviour of mean 
field theories is always the same. The most serious fault of mean field theories lies in 
the neglect of long-range fluctuations of the order parameter. The importance of this 
omission depends very much on the dimensionality of the problem. In problems involv- 
ing one and two-dimensional systems, the results predicted by mean field theory are 
often quantitatively wrong, although they do give a qualitative idea of the behaviour. 



5.2 The triangular lattice model 

First we will look at the development of mean field theory on the triangular lattice, 
which we will use as a basis for our development on the Union Jack lattice. This is 
a well studied area, with many textbooks on the subject, including |LL80j . |KIUH65] 
(summarised in |PB94j ) and |Bax89j . 

Let us consider the triangular lattice Ising model in an external magnetic field, B. 
The Hamiltonian of such a system is given by 

N N 

H = — J ^ cTijai+ij + (Tij(Tjj+i + ajj-aj+ij+i — B ^ cTjj. (5.1) 

Now we focus on one site and consider its expectation value, which at a certain tem- 
perature has value m, 

{cy^) = m (5.2) 

for all i. This m is referred to as order parameter of the system. For the next step, we 
will focus on a particular spin of the lattice ctq. As in our final model all spins will be 
interacting with the parameter m. We can first develop the method for one site with j 
nearest neighbours and then extend it to the entire lattice. For this we take the local 
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version of the Hamiltonian in (5.1), which is 



h (ctq) = -ctq yJ ^ I 

= -ao {qJm + B) - J^q ^ {aj 



m] 



(5.3) 



where q is half the number of nearest neighbours of site 0. This definition of q is such 
that the duphcation of counting of nearest neighbour spins can be ehminated. If we 
disregard the second term in this equation we obtain a non-interacting system. In a 
non-interacting system each spin is an effective magnetic field composed of the applied 
field and an average exchange field due to the neighbours. This means that it will work 
in any number of dimensions, since the lattice configuration only effects the value of q. 
For example q = 2 for a square lattice, q = 3 for both the triangular and simple cubic 
lattice. The magnetisation has to be determined self-consistently from the condition 

m= (ao) = (a,), j = l,2,...N. 



Extending the local Hamiltonian for an non-interacting system to all sites in the matrix 
we obtain the following equation. 



N 



N 



H = —Jqni Gi — B cTj. 



1=1 



i=l 



As before, the partition function can be found from the Hamiltonian 



E 



exp 



N 



{(3Jqm + (3B) Gi 



2 cosh (/3 Jgm + ^B) Zn- 



where Zjq_i 



X^exp 



iPJqm + PB) Oi 



i=l 



and a' denotes a configuration over — 1 sites, (3 = l/{kBT). 



(5.4) 



(5.5) 
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Now we return to looking at the specific case for ctq 



m — 



AT 



1=1 



= Z^^2 sinli {/3Jqm + /3B) ^ exp 

a' 

= Z-^ sinh Jqm + PB)Zn-i 
= t&nh J qm + /3B) . 



{PJqm + PB)J2^^ 

{/3Jqm + (3B) J2 ^^ 



N-l 



i=l 



(5.6) 



To find m, the above equation must be solved numerically for B and T. This equation 
is invariant under the map B — > —B, m — > —m 



— m 



m 



= tanh {-/3Jqm + -f3B) 
= - tanh{-{P J qm + PB)) 
= t&nh {/3 J qm + /3B) . 



(5.7) 



It can be seen from this that for a nonzero external field there is at least one solution 
and sometimes three. For a system with no external field there is always the solution 
m — 0, and two further solutions may exist at ±mo where mo signifies the spontaneous 
magnetisation. As T — >■ 0, tanh (^gJm) — >■ ±1 for m 7^ and mo — >■ ±1. As T 
approaches the critical temperature Tc from T < Tc, |mo(T)| decreases and we can 
obtain its asymptotic dependence from a low order Taylor series expansion of the 
hyperbolic tangent. That is, 



mo = /3g Jmo - - {(3qjf + • 



or 



mo 



(T) = ±V3 



T 



3/2 



T 



1/2 



(5.8) 



From this we can see that the order parameter m approaches zero in a singular fashion 
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as T approaches Tc, vanishing asymptotically as 



mo(T) oc 



T 
T 



1/2 



(5.9) 



It should be noted at this stage that when we compare this with (3.56) that the power 



is \ rather than ^, so will only give a qualitative result, as will be seen later in Figure 



5.3 Extending to the Union Jack 

As we have stated before in Chapter |4| the Union Jack lattice can be thought of as 
two square lattices superimposed on each other. This means that we can develop the 
mean field approximation for this lattice in one of two ways, both of which will be 
shown in this section. First, in the most simple development, we can consider these 
sublattices independently of each other and thus develop a set of partially uncoupled 
equations. In a slightly more complex development, we can approximate the complete 
lattice with a set of coupled equations. That is, equations where there is a dependence 
on both sublattices. In the next section we will perform some analysis of the quality 
of the approximation of both systems. 

5.3.1 Partially uncoupled equations 

The partially uncoupled equations can be quickly developed from the equation for a 
triangular lattice, that is the sublattices have sites with the same number of nearest 
neighbours. As in Chapter |4] we denote the sites with four nearest neighbours as 
and the sites with eight nearest neighbours as cTj. We now simplify the system further 
by taking the interaction strengths along the vertical and horizontal bonds as J, and 
along the diagonals as K. We also include an applied external magnetic field of B. 
Using the above theory for triangular lattices, we can construct equations for each type 
of site: 



{a) = tanh[/3 {2K {a) + B)], 
(r) = tanh[/3(4J((T) +5)]. 



(5.10) 
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In these equations we can see that the number of nearest neighbours is perhaps not 
what we would expect from the previous section. In the equation for (a), we are only 
concerned with interactions between one a and another a, so to avoid duplication we 
have to divide by two. Meanwhile in the equation for (r), we are concerned with 
the number of interactions between a r and a cr, so we do not need to worry about 
duplication as there are no r — r interactions. From this point forward we will re- 



fer to equations (5.10) as the partially uncoupled equations, since both equations are 



dependent on (cr) and so only the equation for (cr) is uncoupled. 



5.3.2 Coupled equations 

While the partially uncoupled equations can be useful, to make the equations uncoupled 
we purposefully ignored the interactions with the (r) in the equation for (cr). Now we 
will develop a system of equations where we consider all the interactions on the a, 
which will be referred to as the coupled equations. Considering the lattice as before, 
we can rewrite the Hamiltonian for the system as 

Again we focus on the expectation values of the sites and put them into our Hamiltonian 
to get 



\ k j / j j k 



The partition function is then 

Z^= J2 /3(2K(f^)+2J(r)+5)^a,+/3(2J((T) + 5)^rfc 

crj,Tk=±l L j k , 

The partition function is defined as the sum over all possible values of cxj and r^, and 
there are N = L"^ sites in total. Let us pick two such spins, say ap and tq, and perform 
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that part of the sum involving these sites. 



^ exp [13 {2K {a) + 2 J (r) + B) ap] 

ap=±l 



= exp [13 {2K {a) + 2 J (r) + B)] 

+ exp [-/? {2K (a) + 2 J (r) + B)] 

= 2cosh[/3(2/s:(a) + 2J(r) + 



Likewise 

exp [p (2 J (a) + B) tq] = 2 cosh [/3 (2 J (a) + B)] . 

rQ=±l 

If Ziv-2 is the partition function for the same system but with N — 2 lattice sites, then 



Zjv = 4 cosh [p {2K {a) + 2 J (r) + 5)] cosh [(3 (2 J (a) + B)] Zn-2- 



Now we compute the magnetisation for ap and tq: 



{ap) = Zjv^ ^ cTpexp 
Now we use the fact that 



/3 (2i^ {a) + 2 J (r) + J]] a^- + ^ (2 J (cr) + ^ 
13 {2K (a) + 2 J (r) + S) J] a,- + ^ (2 J (a) + E) J] 



J2 apexp[p{2K {a) + 2J {t) + B)ap] = 2smh[P {2K (a) + 2J {t) + B)] , 



ap=±l 



J2 rQexp[l3{2J{a) + B)TQ] = 2 sinh [/3 (2 J (a) + 5)] 

TQ=±1 
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to deduce 



(3 {2K {a) + 2 J (r) + B)J2^J + P (2^ (^) + B)J2 



^ ap exp 

0-j,Tfe=±l 

^ sinh [/3 (2fs: (a) + 2 J (r) + fi)] cosh [(3 (2 J (a) + B)] Zn-2 



tanh [f3 {2K {a) + 2 J (r) + B)] , 

J2 rgexp l3{2K{a) + 2J{r) + B)Y,^, + f5{2J{a) + B)Y,rk 

0-j,Tfc=±l L j k . 

^ sinh [/3 (2 J (a) + 5)] Z,v-2 cosh [/3 (2i^ (a) + 2 J (r) + 5)] 
tanh [(3 (2 J (ct) + B)] . 



Since we have (crp) = (a) for all P and (tq) = (r) for all Q we need to analyse the 
equations 



{a) = tanh[(3{2K {a) +2J {t) + B)] 
(r) = tanh [/3 (2 J (a) + 5)] . 



(5.11) 



5.3.3 Wu and Lin adapted 

To test our mean field approximations against the known results of Wu and Lin |WL87t 
IWL89j we will have to rewrite the conditions that they use to classify the type of system. 
To quickly review, for a ferromagnetic system the following conditions hold: 

E, < E2, Es, E^; (5.12) 

for antiferromagnetic systems: 

E2<EuE3,E,; (5.13) 

and for the metamagnetic system: 



Es < Ei,E2,E4 
or ^4 < Ei,E2,E3; 



(5.14) 



5.3 Extending to the Union Jack 



65 



where 



El = J + f + \Ji + J2 + J3 + Jil , 
-E2 = -J - J' + \Ji- J2 + J3- Ji\ 



-E3 = -J + J' + \Ji- J2- J3 + J4 
-Ei = J - J' + \Ji + J2 - J^- Ja\. (5.15) 



Now we rewrite (5.15) with the horizontal and vertical interactions being set to J and 



the diagonal interactions set to K, as before. This produces the following equations: 

-El = K + K +\J + J + J + J\=2K + 4:\J\, 

-E2 = -K - K +\J ~ J + J ~ J\ = -2K, 

-E3 = -K + K+\J-J-J + J\ = 0, 

-Ei = K-K+\J+J-J-J\ = 0. (5.16) 

Using these equations we can rewrite our conditions for the phase of system as: 
to be ferromagnetic: 

K+\J\>0; (5.17) 

to be antiferromagnetic: 

K +\J\ <0 so K < -\J\; (5.18) 
and finally to be metamagnetic 

K + 2\J\ <0 and K > 0. (5.19) 

However \J\ is always positive, so the conditions for metamagnetic systems cannot be 
satisfied. This is due to our simplification to two variables, and means that mean field 
theory will not be able to model metamagnetic systems. At K = it can be seen that 
the Union Jack lattice reduces to the square lattice. 
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5.4 Analysis of the magnetisation 



Having now developed the mean field approximation for the Union Jack lattice, we will 



now analyse the equations (5.10) and (5.11). We will look at the two types of equations 



separately, looking at such quantities as the predicted critical temperatures. Then we 
will go on to look at the graphs of the functions at a constant temperature. As the 
equations for (r) rely on (cr) we will not be discussing them in this section, as they will 
follow those of (cr). Equally as the partially uncoupled equations are adapted forms of 
those from the square lattice, the analysis will be similar and so not included here (see 
|Bax89] ). 



5.4.1 Partially uncoupled equations 



As the equation for (a) in (5.10) only depends on itself no further rearrangement is 
required. For this analysis, and simplicity, we consider a system with no external 
magnetic fields applied, that is 5 = 0. For any temperature T, the corresponding 
value of (a) is determined by the point of intersection of the curves y = x and y = 
tanh [4PK (a)]. The slope of the latter curve varies from the initial value of 4PK at 
the origin and asymptotically approaches the final value of zero. This can be seen 



graphically in Figure 5.1 below. It can be seen that there exists a critical temperature 
above which will have no spontaneous magnetisation. That is the magnetisation will 
be zero. It follows that an intersection of the two curves (other than at the origin) can 
only eventuate if 

A(3K > 1 or T <Tc = (5.20) 

Kb 

where we call the critical temperature. At temperatures below the model acquires 
a non-zero spontaneous magnetisation but at temperature T > Tc no spontaneous 
magnetisation is possible. 



5.4.2 Coupled equations 

Now we move to the coupled equations and perform a similar analysis. However, we 
first need to rearrange the equations to make the equation for a dependent only on the 
value of a. We can do this easily as the equation for r only depends on a so we can 
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Figure 5.1: Graph of the partially uncoupled mean field prediction on the Union Jack 
lattice for a (from eqn. (5.10)), shown in red, against the graph of y = x, shown in 
blue. The intersection of the two lines shows the value of {a) at temperature T. Here 
T = Tc and so the graphs intersect at exactly zero. 



simply substitute it in as shown below: 



(a) = tanh [13 {2K (a) + 2J {t) + B)] , 
(r) = tanh [/3 (2 J (a) + 5)] , 
thus (a) = tanh[(3{2K {a) +2Jtanh [13 {2J {a) + B)]+ B)]. (5.21) 

While this substitution has reduced the number of variables, it has also made the 
analysis more complicated. However at the origin the equation still behaves like the 
linearised version. We can find a critical temperature equivalently by solving 

2/3is: + 4/3V2 > 1, 



or equivalently by solving 

+ 2(3K -1 = 0. 



(5.22) 
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for 13. We put the coefficients into the solution for a quadratic equation 

-2K ±V4J0TTKP 



(3 



8J2 



4J2 

The critical temperature can be found with the following formula 



(5.23) 



4/2 

Tc = , , (5.24) 

ke {-K ± ^K^ + 4 J2) 

This formula produces two roots, and we need to perform some further analysis to find 



which root we should take. We can rewrite (5.23) as a Taylor series expansion for small 
IJI where the ffist term will be 



-K±\K\ (l + S) 



/3 jf, (5.25) 



So for > we ffist take the positive sign, getting 



-K + K (l + Iff) 1 



and if we take the negative sign we get 



B ^ ^ . (5.27) 

^ 4J2 2.P 2K ^ ' 



Now we look when < 0, and again ffist take the positive sign: 



-A'-A'(l + S) -K 



and finally taking the negative sign: 



5.5 Numerical results 
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For ferromagnetic systems, when K > — | J|, we require a positive critical temperature. 



So we take the positive sign, for > using (5.26) and for < using (5.28) 



However, as J approaches zero, the a sublattice becomes a square lattice with the 
critical temperature defined by = 2K/kB- For antiferromagnetic systems we require 
a zero or negative temperature, so we can see for consistency as J approaches zero, we 



must take the negative sign for K < —\J\ (5.29). Here K will always be negative, so 
we only use one equation. 

Again we plot the graph of this function at the critical temperature, shown below 



in Figure 5.2 We can see at the critical point of the system the function follows the 





Figure 5.2: The graph of the coupled mean field prediction on the Union Jack lattice 

X shown in blue. 
Again the 



for a (from eqn. (5.21)) is shown in red against the graph of y 



For this prediction we set the temperature of the system to be again Tc 
intersection of the two lines shows the final value of <t, here at zero. 



line near the origin. As we have seen before with the sites with four bonds, as the 
system passes the critical temperature it falls below this line. 



5.5 Numerical results 

Having developed our mean field approximations we can test them against our known 
results. In this section we will qualitatively compare results from a numerical simulation 
of the mean field approximation (program code given in Appendix [C]) to the known 
results of Wu and Lin |WL87j . While we do not expect the quantitative correlation 
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to be very close, qualitatively the approximation should bring a general idea of the 
character of the magnetisation. We will start by looking at the ferromagnetic square 
lattice and then we will move on to looking at a similar system on the Union Jack 
lattice. 



5.5.1 Triangular ferromagnetic lattice 



First we will look at a triangular lattice where the interactions have strength one 
hundred times the Boltzmann constant and there is no external magnetic field. In 



Figure 5.3 our data points are plotted against the known result. As we can see, the 
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Figure 5.3: Graphs for the triangular lattice Ising model for a ferromagnetic system 
with interaction strength J = 100/cb. Here the solid line shows the prediction of 



Stephenson for nearest neighbour interactions (eqn. (3.57)). The data points show 



the numerical results from the mean field simulation of the system. Note that the 
critical temperature for the mean field curve is higher than the prediction temperature. 
However qualitatively the curves are similar. 



mean field result is approximately the same shape, and shows a phase transition along 
with a critical temperature, of 300 K. That is a critical temperature higher than the 
known result. However the mean field result starts to move away from the positive 
ground state earlier than the analytical result and so has a slightly less steep curve. 
As a qualitative approximation it is not bad, though quantitatively it is quite far from 
the mark. 



5.5 Numerical results 
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5.5.2 Union Jack ferromagnetic lattice 



Next we move on to the Union Jack lattice. In general we have seen a good correlation 
between the theories of Stephenson |Ste70] . Wu and Lin |WL87j . and our intuitive 



expectation in Section 4.3 So as our first example of the results we get from our mean 
field approximation, we will look at the isotropic ferromagnetic case with interaction 
strengths of one hundred times the Boltzmann constant. The graph of the coupled 



equations versus the Wu and Lin's result is shown in Figure 5.4 below. Again we have 
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Figure 5.4: Graphs for the isotropic ferromagnetic system on the Union Jack lattice 
with interaction strengths J = K = lOO/cs. The theoretical predictions of Wu and 
Lin are shown by the solid Hues and the numerical results from the coupled mean field 
equations are shown by the points. Note that the curves are qualitatively the same 
though quantitatively the critical temperature of the mean field results is lower than 
the predictions. 



a good qualitative correlation between the two sets of results, although quantitatively 



the approximation is slightly under the known result. When we compare to Figure 5.3 
we see that it is of about the same quality as the triangular lattice predictions. 



When we look at the partially uncoupled equations we obtain Figure [5~5] below. As 
we can see the correlations of this mean field theory approximation with the known 
results is quantitatively not good. Equally the values for r are larger than those of the 
0", which is in the opposite order to the known result. It is a very rough prediction and 
even when looking at it qualitatively we can see that the curves of the magnetism for 
the two sub lattices are different. In general it would be of more use to use the coupled 
equations as those results are closer to the known theory, so for the rest of the section 
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Figure 5.5: Graphs for the isotropic ferromagnetic system on the Union Jack lattice 
with interaction strengths J = K = lOOA;^. Again the theoretical predictions of Wu 
and Lin are shown by the solid lines. The data points show the numerical results from 
the mean field partially uncoupled equations. Note that the order of the curves for 
the mean field predictions is opposite to the theoretical predictions. Also the critical 
temperature of the mean field is much lower than the theoretical prediction. 



we will concentrate on those equations in (5.11). 

Here we see that the results from the couple equations are more accurate than those 
from the uncoupled equations on isotropic systems. Next we shall look a anisotropic 
system. Of course this allows us to look at a system with a re-entrant phase transition. 
So we will look at the system with horizontal and vertical interactions of J„ = lOOfc^ 
and diagonal interactions oi J = J' = 92kB- The graph we obtain is shown in Figure 



5.6 The result from the mean field has a poor correlation with the prediction of Wu 



and Lin. Qualitatively it can seen that while the sublattice magnetisations are ordered 
in the same way as the predictions, but do not see the re-entrant phase transition at all. 
Interestingly the mean field result here quantitatively is above the critical temperature. 
Also we can see in the graph for the (a) that there is a deviation from a smooth curve. 



5.6 Summary 

In this section, we have shown that the mean field approximation is a powerful tool 
for qualitative prediction of a given system in some instances. It allows us to study 
systems outside the range of the exact solutions as it is not limited by dimensionality. 
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Figure 5.6: Graphs for the anisotropic ferromagnetic system on the Union Jack lattice 
with interactions J = lOOA;^ and K = —92kB- The data points are calculated using 
the coupled mean field equations and the solid lines are the Wu and Lin theorectical 
predictions. Note that the mean field curves do not see the second range of non-zero 
magnetisation. Equally the temperature is higher than all the critical temperatures. 



This allows us to study systems with for example an external magnetic field. It does 
have its limitation however as it is not able to model metamagnetic systems due to the 
simplification used. Also we have seen that to model the Union Jack lattice models we 
require a coupled set of equations to produce results that more closely follow those of the 
known predictions. In addition we see that when the diagonal interactions are negative 



our mean field results do not closely follow the predictions. In the exact solution (4.8) 



the overall power is 1/8, where in the mean field equations (5.11) the overall power is 
1/2. This limits the accuracy of the approximation. The approximation is an iterative 
process, and this could lead to any error being compounded. Thus a defect might lead 
to an instability, resulting in the accuracy of the approximation being poor. 



Chapter 6 

Sampling Methods 



As we have seen from the calculations of Section 3^ and Chapter |4| it is not a practical 
possibility to calculate every state of the Ising model on a lattice of a reasonable size 
as it is NP complete. While one can use an approximation method like the mean field 
theory approach seen in the previous Chapter |5} these generally only have a qualitative 
use and we would like a method that could produce reasonably accurate qualitative 
results. With this in mind in this chapter we will briefly discuss the area of Markov 
Processes |GS92j and how these can be sampled using Monte Carlo Methods to produce 
our results in a reasonable timeframe. While this is going to be a short presentation, 
the reader is directed to books such as |LB05t IBB95j for a more in depth development 
of this theory. Our development however will follow the flow of the development in 
Section 7 of |PB94j . In Chapter [sjwe will use these methods to preform some numerical 
experiments on various Ising model systems. 



6.1 Monte Carlo methods 

In a Monte Carlo simulation |MU49j one does not attempt to simulate the dynamics 
of the system; instead the idea is to generate states ... by a stochastic process 
such that the probability vTj of state i is that given by an approximate distribution (in 
our case the canonical distribution). In a production run of a simulation, states 
are generated and the desired quantity Xi (energy, magnetisation, pressure, etc.) is 
calculated for each state. If the probabilities are correct, then 

i 
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In our case we will calculate the canonical ensemble, so the probability will be given 
by 

TTi = ^e"^^^ where Z = ^e~^^\ 

i 

In the general case a Monte Carlo simulation consists of the following steps: 

1. Choose an initial condition 

2. Select a move 

3. Accept or reject the move using a criterion based on a detailed balance 

4. Repeat steps (2) and (3) until enough data is collected. 

Typically the data from the early part of a run is discarded since the system will not 
have had enough time to reach equilibrium. 

From this definition we now have two questions: 

1. How does the computer generate the states? 

2. How can we make sure the probabilities are correct? 

6.2 Markov processes 

Now let us consider our Ising model. If it has positions then it will have 2^ mi- 



crostates. Using the Hamiltonians given in (3.25), (3.59) and (4.1 ), it is easy to calculate 



the energy of one of these states. If we were to calculate the expected value with the 
canonical ensemble, we would use a random number generator to assign the values of 
spin to each position, weight the contribution of that microstate by e"''^* and repeat 
until the expectation value had converged. This is inefficient since all the states would 
appear with equal probability, including those with such small weight that they do not 
contribute to the thermodynamical average. 

Instead to make the process more efficient, we will be sampling the transition be- 
tween the states. This is because it will be more likely to occur between those of more 
dominance to the system. The process we shall use to generate these states will be a 
Markov Process. 



6.3 The Metropolis-Hastings algorithm 
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If the model starts in a given microstate i, it will move to state j with transition 
probability Pj^i that does not depend on the previous history of the model. If we 
assume the model to be under some fairly general conditions, processes after the passage 
of a transient would produce states with a unique steady-state probability distribution. 
This steady-state probability ttj is an eigenvector, with eigenvalue one, of the transition 
matrix: 

7rj = ^Pj^iTTi. (6.2) 

i 

The steady-state probabilities are unique if the matrix Pj^^i is regular, which means 
that for some integer n all elements of [Pj^^iY are positive and non-zero. Physically, 
this restriction implies that it is always possible to go from one state to any other 
state in a finite number of steps. Exceptions are matrices that are block diagonal, for 
example 





-Pl^2 








-P2-(-l 


-^2^2 
















-P3-(-4 










-P4-(-4 



Since there is no way of going from states 1 or 2 to 3 or 4, the stationary probabihty 
distribution will depend on whether one started with one of the first two states or one 
of the last two. 



6.3 The Metropolis-Hastings algorithm 

Suppose we wish to determine the transition matrix Pi<^j for our Ising model so that 
the steady-state distribution is 

exp -pEi 

7r(z) = (6.3) 

where j3 = \/ {ks T) and Z is the partition function. A possible method of generating 
a sequence of states from an initial state is to pick a site a randomly and attempt to 
flip (change the sign of) its spin. The resulting state (which may be the same state i 
if the attempt to flip a a fails) we call j. Let be the transition probability from i 
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to j. After n steps the transition probability Pf^i{n) is given by 

il,i2,...,in-l 

After many steps the system will approach a limiting distribution 

'^ = exp[-/3{E{m)-E{j)}] 

for all pairs of states mj. We now also require the transition probabilities to be 
normalised 

E^.-- = l (6-4) 



j 



and to obey 



Pj^m _ 7r(j) 



exp[-{Eij)-E{m)}]. (6.5) 



We find that 



Pm^j 7r(m) 

7r(m) = ^P„^j7r(j). (6.6) 

j 

This relation holds, by definition, for any Markov process. The first step above follows 



from the normalisation in (6.4), while the second step involves substituting (6.5). From 



(6.1) we see that (6.6) implies that 7r(m) is a stationary probability distribution of the 



process. Equation (6.5) is called the principle of detailed balance. It can be shown that 
it is a sufficient condition for arriving at the correct limiting probability distribution, 
provided that the process for selecting moves does no contain any traps. That is, it 
should always be possible to get from any given microstate to any other microstate. 

The simplest and most frequently used method of achieving a detailed balance is 
the Metropolis algorithm |MRR"'"53] : 

1. Pick a site a randomly 

2. Compute the energy change AE = E{f) — E{i) that would occur if the spin at 
site a was flipped 

3. If AE < 0, flip the site at site a; if AE > flip this site with probability 
exp(-/3AE) 

4. Repeat steps (1) to (3) until enough data has been collected. 



6.3 The Metropolis-Hastings algorithm 
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An alternative to (3) that is sometimes used is to flip the spin cTq, with probability 



[exp{f3AE) + 1]'^ regardless of the sign of AE. It is easy to see that (6.4) is satisfied in 
both cases for all possible pairs of states. The allowed states therefore will occur with 
the correct frequency if the simulation is run long enough to reach the steady state. 



Chapter 7 

Calibration 



In this chapter we will be calibrating the program written for the project (code given 
in Appendix [B]) by running simulations on a triangular lattice and comparing them to 
the known results given in |Ste64j and |Bax75] . This is an important step as it will 
allow us to see how accurately the program is able to simulate a system quantitatively, 
rather than just qualitatively as seen with mean field theory in Chapter |5| 



7.1 Isotropic system 
7.1.1 Ferromagnetic 

In this first simulation we will investigate an isotropic system with all nearest neighbour 
interactions as one hundred times the Boltzmann constant and a three site interaction 
of zero. This is the simplest system on a triangular lattice, and contains a phase 
transition. This will give us a rough guide as to whether we can continue with more 



complicated models. The resultant graph is shown in Figure |7.1[ We can see from 
this graph that our simulation shows a good correlation with the known results that 
we described earlier. Our simulation tracks the phase transition at under 400 Kelvin 
quite well. The correlation is good for both the magnetisation data and the three- 
site correlator data. There is little noise after the phase transition and the simulation 
results quickly move to a zero average magnetisation. While visible, the noise is of a 
magnitude that is much lower than the rest of the data. 
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Temperature 



Figure 7.1: Graph of numerical simulation results plotted against the theoretical 
predictions of Stephenson |Ste64j and Baxter |Bax75j for a isotropic ferromagnetic 
system on the triangular lattice with interaction strengths Ji = lOOA;^, J2 = lOOfc^, 
J = looks- The theoretical predictions are shown with solid lines. The numerical 
simulations are shown with points of the colour corresponding to the relative prediction. 



7.1.2 Antiferromagnetic 



Having looked at an isotropic ferromagnetic system, the next step is to look at an 
isotropic antiferromagnetic system. According to the conditions set out in the papers 
|Ste64l lBax75j . we should expect results of zero average magnetisations. The graph we 



obtain from the simulations is shown in Figure |7.2[ From the graph we see that the 
simulation data shows a good correlation with the predicted zero magnetisation, and 
appears to have little noise. We also note that there is no evidence of a phase transition 
across our temperature range. This is however a frustrated system. That is, not all of 
the interactions can be in their preferred states, one interaction must always be in a 
parallel state. To further investigate antiferromagnetic systems we have to move on to 
looking at anisotropic systems. 



7.2 Anisotropic systems 
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Figure 7.2: Graph of numerical simulation results plotted against the theoretical 
predictions of Stephenson and Baxter for the isotropic antiferromagnetic system on the 
triangular lattice with interaction strengths Ji = — lOO/cs, J2 = —lOOks, J = — lOOfc^- 



7.2 Anisotropic systems 
7.2.1 Ferromagnetic 

Our investigation of anisotropic systems will start by trying to match the results given 
in Baxter's paper |Bax75j of the three-site interactions on the square lattice. We do 
this by setting the diagonal interaction strength to zero, while choosing values for the 
other two interactions. In our first simulation of these systems we will look at the 
ferromagnetic system where the horizontal and vertical interactions are one hundred 



times the Boltzmann constant. The graph we obtain is show in Figure 7^: Once again 
our simulation has a good correlation with the know results. There is noise around 
the critical temperature, though the results quickly return to the expected value. The 
phase transition appears to be consistent with the known results. 



7.2.2 Antiferromagnetic 

Now we move to a more general anisotropic system in this set of systems, an anti- 
ferromagnetic system where the vertical interactions are of equal magnitude as the 
horizontal interactions but negative. Again we expect the results to show a zero aver- 
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Figure 7.3: Graph of numerical simulation results plotted against Stephenson and 
Baxter theoretical predictions on the isotropic ferromagnetic system on the square 
lattice with interactions Ji = lOOfc^, J2 = IOO/cb, J = 0. 



age magnetisation across the entire temperature range and no phase transitions. The 



graph of this simulation is show in Figure 7_A As expected we can see that the sim- 
ulation data is approximately zero average magnetisation and there is no sign of a 
phase transition in the temperature range. The noise seen is the data is low and seems 
constant across the data set. In this system we do not have the issue of frustrated 
bonds as the lack of the diagonal interactions allows all the interactions to be in their 
preferred state. 

As an extension to the last simulation, we will now perform a simulation of a 
system of an antiferromagnetic system where the vertical and horizontal interactions 
are negative one hundred times the Boltzmann constant and the diagonal interactions 
are positive one hundred times the Boltzmann constant. The graph of this simulation 



is shown in Figure |7.5[ Once again, as with the previous antiferromagnetic systems, 
we see that the data is around the zero average magnetisation. This shows a strong 
correlation to the known results for the triangular lattice, and the noise seems to be 
at a low level. There appears to be no phase transition across the data range, which 
is as expected. In this simulation the interactions will be in their preferred state. 
That is, along the horizontal and vertical bonds the spins will be in an anti-parallel 
arrangement, due to the ferromagnetic coupling of the diagonal bond. 



7.2 Anisotropic systems 
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Figure 7.4: Graph of numerical simulation results plotted against the theoretical 
predictions from Stephenson and Baxter for the isotropic antiferromagnetic system on 
the square lattice with interactions Ji = lOOks, J2 = —lOOks, J = 0. 
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Figure 7.5: Graph of numerical simulation results plotted against the theoretical 
predictions from Stephenson and Baxter for the isotropic antiferromagnetic system on 
the square lattice with interactions Ji = — lOOfcs, J2 = —lOOks, J = WOks- 
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7.3 Overall result 

In the simulations carried out during this chapter we have seen that the data produced 
shows a high correlation with the known results of Baxter and Stephenson. In produc- 
ing accurate results for these systems, confidence can now be taken in future results as 
to their accuracy of modelling the system. In general we can say that the calibration of 
the program has gone well and a useful tool has been developed for the next chapter. 



Chapter 8 

Numerical Simulations 



Having calibrated our program against the triangular lattice, we move on to look at 
the Union Jack lattice. As we saw previously in Chapter [4], the theoretical predictions 
from Wu and Lin |WL89] sometimes give results that intuitively we do not expect. In 
this chapter we will perform numerical simulations on various systems and compare 
the results against Wu and Lin's predicted results, and the critical temperatures again 
against those from Vaks et al. |VLU66j . As in the triangular lattice simulations, we 
will be using a lattice of 100 sites by 100 sites in our program. The simulation results 
obtained will be for a finite lattice, while the theoretical predictions are for an infinite 
lattice. This lattice size is chosen as it is small enough to have a reasonable run time, 
while being large enough to suppress the finite size effects. 



8.1 Isotropic ferromagnetic 



Following the structure of our analysis in Section 4^, we will start our numerical 
simulations with the isotropic ferromagnetic system. Here the interactions for our 
system will be J = J' = J„ = lOOks- When we plot our simulation results against 



the predicted results of Wu and Lin we obtain Figure |8.1 below. As we can see 



our simulation data has a high correlation with the prediction functions. As we have 
discussed before, Wu and Lin's result follows that of Vaks et al, and so we can say 
that our simulation data also agrees with their phase predictions. There is some noise 
around the critical temperature, though is of a small magnitude when compared to the 
other results. After the noise part of the data, the points again follow the prediction 
with a high correlation. 
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Figure 8.1: Graph of numerical simulation results plotted against theoretical predic- 
tions of Wu and Lin for an isotropic ferromagnetic system on the Union Jack lattice 
with interactions Ji = WOks, J2 = WOks, J3 = WOks, Ja = 100/cbi J = lOOfcs, 
J' = lOOfcs. Here the theoretical predictions are shown with the solid lines and the 
numerical results are shown with the points. 

8.2 Anisotropic metamagnetic 

We move on to the first system in which we saw a disagreement between the results 
of Wu and Lin and the configuration we expect. Our next simulation will be on an 
anisotropic metamagnetic system where J„ = lOA;^, J = lOOA;^ and J' = lOOks- The 
graph we obtain when we plot our simulation results against the predictions of Wu and 



Lin is shown in Figure [872| below. Our simulation results show a poor correlation to the 
prediction functions, with the average spin across the temperature being zero. From 



our analysis in Section |4.3| we concluded that the system should be antiferromagnetic 
when there is no external magnetic field. The results have a high correlation to an 
antiferromagnetic system, and so back up this analysis. 



8.3 Anisotropic antiferromagnetic 

Next we look at the anisotropic antiferromagnetic system. The system will have hor- 
izontal and vertical interactions of J„ = lOOks and diagonal interactions of J = J' = 
— looks as before. The results are plotted against Wu and Lin's predictions in Figure 



8.3 below. The correlation between the simulation data and Wu and Lin's prediction is 



8.3 Anisotropic antiferromagnetic 
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Figure 8.2: Graph of numerical simulation results plotted against theoretical predic- 
tions of Wu and Lin for the anisotropic metamagnetic system on the Union Jack lattice 
with interactions J„ = lOks, J = WOks and J' = — lOOfc^. Note that the simulation 
results do not follow the curves of the prediction functions. 
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Figure 8.3: Graph of numerical simulation results plotted against theoretical pre- 
dictions of Wu and Lin for an anisotropic antiferromagnetic system on the Union 
Jack lattice with interactions Ji = lOOA;^, J2 = lOOfc^, J3 = lOO/c^) -A = lOOA;^, 
J = — lOO^B, J' = —looks- Note that the simulation results do not follow the curves 
of the theoretical predictions and are physically possible. 
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again poor. Vaks et al. predict that this system should have average spin of zero. As 
we can see this confirms the results of the simulation. Overall the noise in the results 
is low and does not show any signs of detecting any non-zero magnetisation. 



8.4 Anisotropic ferromagnetic 

So far in simulations we have focussed on systems with only one dominant system. 



In Section 4.3 we saw that we can have an anisotropic ferromagnetic system where 
there is a re-entrant phase transition in the temperature range. The system we studied 
there was one with horizontal and vertical interactions of J„ = lOOfc^ and diagonal 
interactions of J = J' = —92kB. The graph of our simulation results plotted against 



the prediction functions of Wu and Lin is shown in Figure 8.4 We can see that our 
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Figure 8.4: Graph of numerical simulation results plotted against theoretical predic- 
tions of Wu and Lin for an isotropic ferromagnetic system on the Union Jack lattice 
with interactions Ji = lOOfce, J2 = WOks, J3 = WOkB, Ji = lOO/c^, J = —92kB, 
J' = —92kB- Note that the simulation results have good correlation up to the first 
critical temperature, but do not show the re-entrant phase. 



simulation results show good correlation with the predictions of Wu and Lin up to 
the first phase transition. Note that our simulation results do not show the re-entrant 



phase transition. As we discussed in Section |4.3| the re-entrant phase transition is 
between a disordered antiferromagnetic phase and an ordered antiferromagnetic phase. 
Due to the set up of the simulation program we are not able to see the difference 
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between the different antiferromagnetic phases. Compared to the prediction functions 
of Wu and Lin, we can see that our simulation resuhs do not follow their predictions 
for the ordered antiferromagnetic phase, although do show an average magnetisation 
of zero. 



Next in our study of the systems in Section 4^, we will now go on to look at 
the anisotropic ferromagnetic system with horizontal and vertical interactions of J„ = 
— lOO^B and diagonal interactions oi J = J' = 92k b- The graph of this simulation is 



show in Figure 8.5 Here we see that our simulation results have good correlation with 
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Figure 8.5: Graph of numerical simulation results plotted against theoretical pre- 
dictions of Wu and Lin for an anisotropic ferromagnetic system on the Union Jack 
lattice with interactions Ji = — lOO/c^, J2 = — lOOA;^, J3 = — lOOfcs, J4 = —WOks, 
J = 92kB, J' = 92kB- Note that the simulation results shows high correlation with 
both prediction curves. 



the predicted functions of Wu and Lin. Our simulation also shows the slight overall 
magnetisation after 200 Kelvin, up to the critical temperature. There is noise after the 
critical temperature as the two sets of data converge to the average zero spin. We have 
seen this in all our ferromagnetic systems and it is not significant given the range of 
the rest of the data. So we have agreement with Wu and Lin's prediction. 

Following on in our studies, we will look at the anisotropic systems where rotation 
produces different graphs. We first look at the system with horizontal interactions 
of Ji = J3 = lOOfc^/O.Q^, vertical interactions of J2 = = 100^^/0.9 and diagonal 
interactions oi J = J' = IOO/cb, and the similar system rotated through 90 degrees. 
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The graphs of our simulation results are shown in Figure [8^ below. As we see when 
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Figure 8.6: Graphs of numerical simulation results plotted against theoretical predic- 
tions of Wu and Lin for an anisotropic ferromagnetic system on the Union Jack lattice, 
(a) shows the system with interactions Ji = J3 = IOO/cb/0.9^, J2 = Ja = 100^^/0.9 
J = J' = lOO/cfi. (b) shows the rotated system with interactions Ji = J3 = 100A;b/0.9, 
J2 = Ja = 100A;b/0.9^ J = J' = lOOA;^. Note that although the prediction functions 
are different, the simulation results are identical. Also there is poor correlation with 
the T-sublattice prediction. 



we compare the simulations results, the data forms similar curves and has a similar 
phase transition at equal critical temperatures. In comparison to the predictions of 
Wu and Lin, we see that at the higher temperatures the data follows all three curves 
with good correlation. At lower temperatures, below about 200 Kelvin, we see that 
the a prediction still has good correlation for both systems. This confirms that the 
prediction for r (and subsequently the mean magnetisation) is not correct below this 
temperature. It also shows that rotating the lattice should not have an effect on the 
results of the system. 

A further result can be seen if we now take a system similar to the previous example 
but with negative diagonal interactions. For an example of this type of system we will 
look at the system with horizontal interactions of Ji = J3 = 100A;b/0.9^, vertical 
interactions of J2 = Ja = lOOfc^/O.Q and diagonal interactions of J = J' = — 100/cb, 
and the same system rotated through 90 degrees. The graphs of our simulation results 
are shown in Figure 8/7, Again we see as in Figure 8.6| that the simulation results 
are very similar to each other. However, in this case we see that all three simulation 
results have a range of lower results at low temperatures. After these lower results the 
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Figure 8.7: Graphs of numerical simulation results plotted against theoretical predic- 
tions of Wu and Lin for an anisotropic ferromagnetic system on the Union Jack lattice, 
(a) shows the system with interactions Ji = J3 = WOks/O-d^ , J2 = J4 = 100^:^/0.9 
J = J' = —lOOkB- (b) shows the rotated system with interactions Ji = J3 = 
IOO/cb/0.9, J2 = Ji = lOOfcs/O.g^ J = J' = -lOO/cB. Note that although the pre- 
diction functions are different, the simulation results are identical. Also there is slight 
deviation at low temperatures of the simulation results. 

simulations then move up again to the prediction curves, following the a prediction in 
both cases. We do see that in both cases the simulation results and predicted results 
show the same critical temperature and phase transition. This oddity maybe due to 
the simulation being performed on a finite system, while the theoretical predictions are 
for an infinite lattice. 



8.5 Magnetic fields 

One advantage of using a numerical simulation, is that we can study systems in external 
magnetic fields. Metamagnetic systems have the property that when we add a external 
magnetic field we obtain results that show average magnetisation of a higher magnitude 
than that of the applied field. In all the following systems, the base system is one 
where the square lattice interactions are J„ = lOks and the diagonal interactions are 
J = —J' = looks, our standard metamagnetic system. In our first simulation we will 
apply a magnetic field of strength equal to the Boltzmann constant. The graph is show 



in Figure 8.8 With this magnetic field we can see already that at low temperatures that 



there is now a non-zero average magnetic spin. Also the sublattices now have distinct 
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Temperature 



Figure 8.8: Graph of the simulation results for the anisotropic metamagnetic system 
on the Union Jack lattice in an external magnetic field of strength ks- Note that at low 
temperatures there is non-zero magnetisation of greater magnitude than the external 
magnetic field. 



curves at these temperatures, rather than being in the disorder antiferromagnetic state. 
The average magnetisation is small in comparison with the possible values for the spin 
states. So our next step is to increase the magnetic field to a strength ten times 



the Boltzmann constant. The graph we obtain is shown in Figure 8.9 below. We 
can see again that a non-zero average magnetisation has resulted. We see also that 
the a prediction is now showing more of a sign of a ferromagnetic curve, while the 
T prediction is a less sharp curve. Above 200 Kelvin we see a slight increase in the 
magnetisation suggesting that this is the critical temperature where the system is going 



into a disordered phase. When we look back to Figure |45] we see that this is also the 
point where Wu and Lin's prediction hits the x-axis. 

For completeness we will now look at the same metamagnetic system in an external 



magnetic field of lOO/c^- The graph of the results we obtain is shown in Figure 8.10 



Here we see that the results are showing more of a ferromagnetic shape. Again above 
200 Kelvin the sublattice results are converging to a central value, and all the curves 
seem to be related to each other. Below this temperature we see that the r-sublattice 
is showing more of a ferromagnetic phase, although the curve is not sharp. Meanwhile 
the (T-sublattice has a shape suggesting a ferromagnetic phase at low temperatures but 
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Figure 8.9: Graph of the simulation results for the anisotropic metamagnetic system 
on the Union Jack lattice in an external magnetic field of strength lO/cs. 
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Figure 8.10: Graph of the simulation results for the anisotropic metamagnetic system 
on the Union Jack lattice in an external magnetic field of strength lOO/c^- 
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moves up to a maximum at 200 Kelvin. 

8.6 Summary 

We have seen in this chapter that our simulation results help support our findings in 
Chapter |4j We have good correlation between our simulation results in isotropic fer- 
romagnetic and most anisotropic ferromagnetic systems and Wu and Lin's predictions 
|WL87t IWL89j while in a ferromagnetic phase. However in terms of the r-sublattice 
prediction we saw that our simulations support that the results should not change by 
a 90 degree lattice rotation. We also saw that our simulation agrees with Valks et al. 
predictions |VL066j for the antiferromagnetic systems. In addiction we showed that 
the re-entrant phase transition can not be seen from average magnetisation alone. 

In addition we showed that our simulation can be used to further investigate system 
properties that occur in external magnetic fields. This shows that we have a powerful 
tool for numerical simulation of the Union Jack lattice. 



Chapter 9 

Conclusion 



Over the course of the thesis we have seen two main results. The re-entrant phase 
transition can not be seen when looking at the average magnetisation results. In 
addition we have seen that the prediction functions of Wu and Lin in their papers 
|WL87l IWL89j produce some anomalous results. In this section we will review the 
results from the different approaches we undertook and discus the disagreements from 
Wu and Lin's results that were found. At the end of the section we will discuss the 
additional conditions that when applied allow the prediction functions of Wu and Lin 
to show a greater agreement with known theories. 

Our aim was to investigate the properties of the model and specifically the re-entrant 
phase transitions. From our investigations we have seen that these re-entrant phase 
transitions can not be seen when only looking at the average magnetisation results. 
This is due to the transition being from an unordered antiferromagnetic phase to an 
ordered antiferromagnetic phase. While this can been seen in the Vaks et al. paper 
|VL066j . this can be seen from other order parameters, the average magnetisation 
of an unordered antiferromagnetic phase and an ordered antiferromagnetic phase is 
identically zero. This was further confirmed by our numerical simulations, which also 
did not see these re-entrant phase transitions due to being focused on the average 
magnetisation of the system. However we have seen that the predictions of Wu and 
Lin can be used to find the critical temperatures of these phase transitions, and agree 
with the calculations of Vaks et al. 

In addition, we have seen that Wu and Lin's prediction for the a-sublattice given 
in |WL87j requires additional conditions to agree with our numerical simulations. It is 



possible to classify the phases of the system by examining the 7 terms of equation (4.9 ). 
Wu and Lin's predictions under the current conditions produce non-zero magnetisations 
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for non-ferromagnetic systems. However if we impose the conditions that we only use 



the prediction formula (4.8) when 71 < or 71 72 73 74 > and it is zero outside 
those conditions, their results now work for all systems. 

When we look at Wu and Lin's prediction for the r-sublattice, given in |WL89j . 
we see for the general anisotropic lattice there are some additional issues. When the 
interactions on the square lattice are equal, by applying the conditions above we can 
eliminate the physically impossible results that we saw in the antiferromagnetic and 
metamagnetic phases. This is in agreement with the earlier work of Lin and Wang 
|Lin88j . When we consider unequal vertical and horizontal interactions on the square 
lattice we see that physically plausible results are only possible when the horizontal 



interactions are larger than the vertical interactions. From Figure 8^ we see however 
that our simulation results show that there is still some disagreement at low tempera- 
tures. While this suggests that further conditions are required for this prediction, they 
would be more involved than those for the a-sublattice. For the current prediction, 
imposing the condition that Ji = J2 = J3 = Ja, together with the conditions imposed 
on a would bring agreement with the simulation results. 

From our investigation of mean field theory, we saw that it is a poor quantitative 
predictor for the systems. In addition it has the limitation of not being able to look at 
metamagnetic systems. Qualitatively however, it does allows a prediction of the general 
behaviour of symmetric systems without a re-entrant phase transition. Again the mean 
field prediction does not see re-entrant phase transitions in terms of the magnetisation. 
In these systems its prediction of critical temperature is generally in line with the third 
phase transition critical temperature. The results were useful for confirming those the 
results we found in our theoretical analysis. 

In conclusion, in the thesis we have not been able to investigate the re-entrant phase 
transition using the average magnetisation. However, we have developed an accurate 
simulation program for finite systems for both the triangular and Union Jack lattices. 
This simulation program has allowed us to investigate inconsistencies in the theoretical 
predictions of Wu and Lin |WL87l IWL89] . We have also shown that it is possible with 
additional conditions to make these predictions model the systems more accurately. 

Future research could extend this work to look more closely at the rotational vari- 
ance that is observed in the r-sublattice prediction. Intuitively the predictions should 
be invariant under rotation of the lattice, as shown in the simulation results, and this 
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disagreement requires further investigation. Equally as we saw in Figure 8/7 that there 
is an oddity in the results at low temperatures, and this should be investigated to 
determine if it is affected by the lattice size. 

As well as looking to the current lattice configuration, further research could take 
place to see if the results found occur in other configurations. Research into the hexago- 
nal and double hexagonal lattices may be interesting to determine if a re-entrant phase 
is present. In addition in terms of the Union Jack lattice, an investigation into mixed 
spin lattices such as those of Strecka et al. |Str06t ISvDOGj may bring insight into both 
the re-entrant phase transitions and the metamagnetic systems. 

In terms of the simulation program, other Monte Carlo methods could be investi- 
gated to improve its efficiently. After a conversation with Dr Tim Garoni of the Uni- 
versity of Melbourne at the Australian Mathematics Society meeting in 2010, the algo- 
rithm of Swendsen and Kotecky |SW87j , which has been improved by Wang, Swendsen 
and Kotecky in |WSK90j . was suggested to be a more efficient method. As this is a 
cluster method, its use on a Union Jack lattice would need to be investigated. If the 
algorithm could be used on this lattice type, it would overcome the problem of non- 
ergodic data in low temperature antiferromagnetic systems. It would also allow the 
simulation to not be affected by the critical slowing down effects around the second- 
order phase transitions. 
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Appendix A 

Constants and Conversion Factors 



Boltzmann Constant (/cs) 

Avogadro's number , 

1 electron-volt (eV) 

Kelvin 



8.617343 X lO-^cVK"^ 

6.02 X lO^^mol-^ 

1.60 X lO-i'^J 

or 1.1605 X lO^K 
-273.16°C 
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Programming I 



B.l Programmer's guide 

The broad purpose of this code is to simulate the Ising Model on a two-dimensional 
lattice with specific configurations. It does this by first asking the user for the initial 
conditions of the system along with the configuration of the lattice, the temperature 
range required and how many initial points to ignore. The program then simulates the 
model over randomly chosen temperature points, and outputs the data to a text file. 
The program has the ability to be split the output into sublattices if required. 

B.1.1 Flow 

Initially the user is asked for the type of lattice, the interaction strengths, external 
magnetic field strength they want the system to be running with, the number of points 
to ignore and the temperature range they are interested in. After all the conditions 
have been entered the system is initialised with all the sites being set to the "UP" (+1) 
position. The initial energy is then calculated for the system. Then a temperature is 
picked at random within the range given by the user, and the system is simulated at this 
temperature. Initially the first set of data points are ignored up to the level required 
by the user. Then the next points are averaged and stored as the average spin value for 
that temperature. The user is then told the overall average spin. The system is reset to 
the " UP" position for every site and is simulated again at a random temperature, not 
equal to the past temperatures. After 200 data points have been stored in the array, 
it is then outputtcd in a suitable format to a text file. The program then ends, telling 
the user how long the simulation took. 
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Figure B.l: Flowchart of the numerical simulation program 

B.2 Code listings for main program 
B.2.1 isling2dt.h 

1 %auto— ignore 



// standard C headers 
#:include <stdio.h> 
#include <math . h> 
^include <stdlib.h> 
#include <time.h> 

// Some useful system variables 
enum { 



10 


VP= 1 , 


11 


DOWtt -1, 


12 


ITERATIONS= 10000, 


13 


MODELS= 1000, 


14 


(X)ND= 1000, 


15 


length= 100, 


16 


height= 100, 


17 


array_len= 200 


18 }; 





B.2 Code listings for main program 
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19 

20 #define BOLT 8.617343*pow(10 , -5) /*The boltzman contant*/ 
21 



22 


typedef 


struct 


actions { 


23 




double 


jh; 


24 




double 


jh2; 


25 




double 


jv ; 


26 




double 


jv2 ; 


27 




double 


jdu ; 


28 




double 


jdd; 


29 




double 


h; 


30 




double 


ku ; 


31 




double 


kd; 


32 




char type ; 


33 


} ACT; 







34 

35 typedef struct coordinate { 

36 float temp; 

37 double p; 

38 double q ; 

39 double r ; 

40 double tripu ; 

41 double tripd ; 

42 struct coordinate *next; 

43 } COORD; 
44 

45 //Prototypes of our functions 

46 int sctup_chain (char *) ; 

47 double set up _gen ( char *, char *); 

48 void i n i t i a 1 i s e _ c h a i n ( int [ length ][ height ] ) ; 

49 void * safe _malloc ( siz e _t , char *) ; 

50 void *safc_rcalloc (int *, sizc_t , char *); /* Does the same for realloc */ 

51 void metropolis ( int [ length ][ height ] , float, ACT, double [3], int [3]); 

52 double energy (int [ length ][ height ] , ACT); 

53 double energyb(int [ length ][ height ] , ACT, int, int, int); 

54 void wr i t e _2 d _ar r ay (char *, ACT); 

55 void set-temp ( float , float); 

56 void free_cond (void); /* Frees the memory taken by linked list */ 
57 

58 extern COORD *holl ; 
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B.2.2 isling2dt.c 



1 %auto— ignore 

2 /* We link our files with a standard header */ 

3 #include " i s 1 i ng 2 d t . h" 
4 

5 OOORD *holl= NULL; 
6 

7 int main() { 

8 // intialising variables 

9 int **chain; // for our lattice 

10 int i , point , count , cycle , ignore , liours , minutes , seconds , * lattice ; 

11 int t imeleft ; 

12 time_t timer ; 

13 char input , save ; 

14 aX)P!D * coord ; 

15 double j , k ; 

16 float min , max, minmax, percent; 

17 ACT vars ; 

18 double * split; 

19 char filename [ 1 0] ; 

20 min= 0; 

21 max= min ; 

22 // some intial values incase of square lattices 

23 vars . ku= 0; 

24 vars . kd= 0; 

25 vars . jv= 0; 

26 vars . jh= ; 

27 vars . jdu= 0; 

28 vars . jdd= 0; 

29 // Finding out the type of lattice required 

30 while ( vars . type != 'S' &fe vars . type != 'T' &k vars . type != 'U' &&: vars . type ! = 

'N' &fe vars. type != 'H' && vars . type != 'D ' ) { 

31 printf("What type of lattice would you like to study? \n ([S]quare, [T 

] riangular , [U] nion Jack , [N] ext Nearest , [H] exagonal , [D] ouble 
Hexagonal ) " ) ; 

32 fscanf ( stdin , "%ls" , &vars.type); 

33 vars . type= toupper ( vars . type ) ; 

34 } 

35 while(input != 'A' && input != 'I') { 

36 printf("Is the lattice [ I ] somorpic of [A] ntimorphic ? (I or A) "); 

37 fscanf (stdin , "%ls" , feinput ) ; 

38 input= toupper ( input ) ; 

39 } 

40 // Finding out the attraction forces between particles 

41 if (input = 'I ') { 

42 // for the is otopic lattice 

43 j= setup_gen (" the attraction \n between two particle", "real"); 

44 vars.jh= j; // if isotopic all j values are the same 

45 vars . j v= j ; 
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46 vars . jh2= j ; 

47 vars . j v2= j ; 

48 if(vars.type != 'S' &fe vars . type != 'H' &fe vars . type != 'D') { 

49 // specfic qualities that are different from square lattice 

50 vars . jdu= j ; 

51 k= setup_gen (" the attraction between three particles", "real") 

52 vars.ku= k; // as all values of k are the same too 

53 if (vars. type = 'U' || vars . type = 'N') { 

54 // specfic qualities for Union Jack and NN lattices 

55 vars . j dd= j ; 

56 vars . kd= k ; 

57 } 

58 } else if(vars.type = 'D') { 

59 k= setup_gen (" the attraction between two particles second 

hexagon", "real"); 

60 vars . ku= k ; 

61 vars . kd= k ; 

62 } 

63 } else if(input = 'A') { 

64 // for the ant it o pic lattice 

65 vars.jh= setup_gen ( " the attraction \n between two particle , horizontal 

" , " real") ; 

66 vars.jv= setup_gen (" the attraction \n between two particle , vertical", 

" real" ) ; 

67 vars.jh2= vars.jh; 

68 vars.jv2= vars.jv; 

69 if (vars. type != 'S' &fe vars . type != 'H' &fe vars . type != 'D') { 

70 // specfic qualties that are different from square lattices 

71 vars.jdu= setup_gen (" the attraction \n between two particle , 

diagonal up", "real"); 

72 vars.ku= sctup_gcn (" the attraction \n between three particles 

up" , " real" ) ; 

73 if (vars. type = 'U' || vars . type = 'N') { 

74 // specfic qualities for the Union Jack lattice 

75 if (vars. type = 'U') { 

76 printf("Does Jl = J3 AND J2 = J4?: "); 

77 fscanf (stdin ,"%ls",&save) ; 

78 printf ("\n" ) ; 

79 save= toupper ( save ) ; 

80 if(save = 'N') { 

81 vars.jh2= setup_gen (" the attraction \n 

between two particle , horizontal 
(J3)" , " real" ) ; 

82 vars.jv2= setup_gen (" the attraction \n 

between two particle , vertical ( 
J4)" , " real" ) ; 

83 } 

84 } 
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85 vars.jdd= setup_gen (" the attraction \n between two 

particle , diagonal down", "real"); 

86 vars.kd= sctup_gcn (" the attraction \n between three 

particles down", "real"); 

87 } 

88 } else if(vars.type = 'D') { 

89 vars.ku= set up_gcn ( " t he attraction between two particles , 

horizonatal second hexagon" ," real" ) ; 

90 vars.kd= setup_gen (" the attraction between two particles , 

vertical second hexagon", "real"); 

91 } 

92 } 

93 // Finding the overal field 

94 vars . h= setup_gen("the overall field", "real"); 

95 // Setting the number of intial values to ignore 

96 ignore= ( int ) set up_gen ( " how many tens of inital values to ignore", "integer"); 

97 while(ignore >= MODELS || ignore < 0) { 

98 // in case they got it wrong 

99 pr in t f (" Error : too many values ignored. Please enter a positive number 

less than %d\n" , MODELS); 

100 ignore= ( int ) setup_gen ("How many inital values to ignore", "integer"); 

101 } 

102 while (min ^= max) { 

103 //setting the minimum and maximum temperatures 

104 min= ( float ) set up _gcn (" minimum temperature", "real above zero"); 

105 max= ( float ) setup-gen ("maximum temperature" , "real above zero and min" 

); 

106 if (min < 0) min= 0; //As we are working in Kelvin it makes little 

sense to have negative numbers 

107 if (max < 0) max= 0; 

108 } 

109 //to save retyping if the temperatures are inputted the wrong way round 

110 if (min > max) { 

111 minmax= min ; 

112 min= max; 

113 max= minmajc; 

114 } 

115 printf (" Enter output filename: "); 

116 if (fscanf (stdin , "%s" , fefilename ) != NULL) { 

117 printf ("\n" ) ; 

118 fflush (stdin) ; 

119 } 

120 timer= time(NULL); // starting the timer 

121 // seeding the random numbers with the current time 

122 srand ((unsigned ) time (NULL) ) ; 

123 // creating our lattice 

124 chain= safe_malloc ( length * sizeof(int [height]) , "Allocating Memory to Chain") 

125 // The simulation starts here 

126 i n i t i al i s e _ c h a i n ( chain ) ; 
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127 coord= hoU ; 

128 if(vars.type = 'U') { 

129 lattice= safc_iiianoc (5 * s iz eof ( int )," Holding array for spin"); 

130 split= safe_malloc (5 * sizeof (double) ," Holding array for average spin" 

); 

131 } else { 

132 latticc= safc_malloc (3 * sizeof ( int )," Holding array for spin"); 

133 split= safe_malloc (3 * sizeof (double) ," Holding array for average spin" 

); 

134 } 

135 for(point= 0; point < array_len ; point++) { 

136 // getting a temperature 

137 set _temp (/* coord , */ min , max); 

138 // starting our average anew 

139 coord= holl ; 

140 lattice [0]= length*height ; 

141 if ( vars . type = 'T') { 

142 lattice [1]= 2* (((length -1)*( height -1) )+(( length -l)+(height -1)) 

); 

143 } else if (vars. type = 'U') { 

144 lattice [1]= length* height * . 5 ; 

145 lattice [2]= length* height * 0. 5 ; 

146 lattice [3]= ((( length -1) *( height -1) )+(( length -l) + (height -1) ) ) ; 

147 lattice [4]= ((( length -1) *( height -1) ) +((length -l)+(height -1) ) ) ; 

148 

149 } 

150 // setting the chain to "UP" 

151 if (point != 0) i n i t i a 1 i s c _c h a i n ( chain ) ; 

152 for(cycle= 0; cycle < MODELS; cyclc++) { 

153 for ( count =0; count < ITERATIONS; count++) { 

154 //running the metropolis algorithm 

155 metropolis (chain , coord— >tcmp , vars, split, lattice); 

// Performing the algorithm 

156 if (cycle ^= ignore &fe count = 0) { 

157 split[0]= (double) lattice [0] /( length*height) ; 

158 if(vars.type = 'T') { 

159 split [1]= (double) lattice [l]/(2*((( 

length -1)*( height -1) )+(( length -1) 

+ (height-l)))) ; 

160 } else if (vars. type ^= 'U') { 

161 split [1]= (double) lattice [l]/( length* 

height *0.5) ; 

162 split [2]= (double) lattice [2] /( length* 

height *0.5) ; 

163 split [3]= (double) lattice [3]/(2*((( 

length -1)*( height -1) )+(( length -1) 
+(height-l)))) ; 

164 split[4]= (double) lattice [4] / (2* (( ( 

length -1)*( height -1) )+(( length -1) 
+(height-l)))) ; 
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165 } 

166 } 

167 } 

168 percent= (( float )( cycle+l+(MODELS* point ))/( float ) (MODELS* 

array _Ien ) ) ; 

169 seconds= ( time (NULL)— timer ) ; 

170 timclcft= (long int )((( time (NULL)— timer )/ percent )— seconds ) ; 

171 seconds= timeleft%60; 

172 minutes= abs(timeleft /60)%60; 

173 hours= abs ( timeleft /(60*60) ) ; 

174 printf ("%.2f %% done , time left: %.2d hours, %.2d minutes, %.2 

d seconds \r" , percent * 100 , hours, minutes , seconds ) ; // To 
show the program is still working 

175 } 

176 coord->p= (double) split [0]/(double) ((MODELS-ignore)*ITERATIONS) ; // 

storing the average spin 

177 if(vars.type = 'T' || vars . type = 'U') { 

178 coord->q= (double) s plit [ 1 ]/( double) ( (MODELS-ignore ) ITERATIONS 

); 

179 } 

180 if (vars . type = 'U') { 

181 coord->r= ( double ) s p 1 i t [ 2 ]/( double )( (MODELS-ignore ) ITERATIONS 

); 

182 coord->tripu= (double ) s p 1 it [ 3 ]/( double )( (MODELS-ignore ) * 

ITERATIONS) ; 

183 coord->tripd= ( double ) s p 1 i t [4 ]/( double )( (MODELS-ignore ) * 

ITERATIONS) ; 

184 } 

185 printf (" \nThc average is %.2 f \n" , coord— >p) ; // a check to see if it is 

resonable 

186 printf("The temperature is %.2 f \n" , coord— >temp) ; 

187 printf (" Point : %i\n" , point); // where we are up to 

188 } 

189 write_2d_array (filename, vars); // outputing to our data file 

190 free_cond(); 

191 printf ("Now run maple to plot graph\n"); // handy hint to the user 

192 free ( chain ) ; 

193 free ( lattice ) ; 

194 free(split); 

195 seconds= ( time (NULL)-t imer ) %60; 

196 minutes= ( abs (( time (NULL)— timer— seconds ) /60) ) %60; 

197 hours= abs (( time (NULL)-timer-seconds -(minutes *60) )/ (60*60) ) ; 

198 printf (" Runtime : %ld hours, %ld minutes, %ld seconds\n" , hours, minutes, 

seconds); // how long it too 

199 return 0; // finish 



200 } 
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B.2. 3 setup2d.c 

1 %auto— ignore 

2 /* We link our files with a standard header */ 

3 #include " i s 1 i ng 2 d t . h" 
4 

5 void * safe _m alloc (size_t size , char *location) { 

6 /* This function performs a malloc operation but checks to 

7 see if there is enough free memory for this to be written to. */ 

8 void * ptr ; 

9 ptr= malloc { size ) ; 

10 if (ptr = NULL) { 

11 printf("Out of memory at function: %s\n" , location ) ; 

12 exit(-l); 

13 } 

14 return ptr ; 

15 } 
16 

17 void *safe_realloc(int * array , size_t size, char *location) { 

18 /* This is like the safc-malloc but for realloc. */ 

19 void * ptr ; 

20 ptr= r e al 1 o c ( array , size); 

21 if (ptr = NULL) { 

22 printf("Out of memory at function: %s\n" , location ) ; 

23 exit(-l); 

24 } 

25 return ptr ; 

26 } 
27 

28 int setup_chain (char *quant) { 

29 /* This function is a left over from a previous program for 

30 the one dimensional Ising model. A two dimensional version 

31 is being looked into.*/ 

32 int num; 

33 print f (" Enter the %s of chain required: (0,1,...) ", quant); 

34 fscanf (stdin , "%d" , &num) ; 

35 printf ("\n" ) ; 

36 return num; 

37 } 
38 

39 void i n i t i a 11 s e -C h a i n ( int chain [ length ][ height ] ) { 

40 /* This function take the chain and set all site in the "UP" 

41 position */ 

42 int i , j ; 

43 for(j=0; j < height; j++){ 

44 for(i=0; i < length; i++){ 

45 if (chain [ i ] [j ] != UP) chain [ i ][ j ]= UP; 

46 } 

47 } 

48 } 
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49 

50 double setup_gen (char *item , char *units) { 



51 /* This is the basic input function , as so many of our variables 

52 have similar questions I decided to write a function for them*/ 

53 double i ; 

54 printf (" Enter the value of the %s : (%s) ", item , units ) ; 

55 if (fscanf (stdin , "%lf", &i ) ) { 

56 printf ("\n") ; 

57 } else { 

58 printf (" \n EEmOR: No input found. \n" ) ; 

59 clearerr(stdin); 

60 fflush ( stdin ) ; 

61 setup_gen ( item , units ) ; 

62 fflush (stdin) ; 

63 } 

64 return i ; 



65 } 

66 

67 void write_2d_array (char output [] , ACT action) 

68 /* Take the arrays in x and y both of length array.len and write a file in 

69 maple output format to the file named filename */ 

70 { 



71 int i ; 

72 int m= array.len ; 

73 COORD * coord; 

74 FILE *fptr ; 

75 coord= holl ; 

76 fptr= fopen (output, "w" ) ; 

77 fprintf (fptr, "p;=[\n"); 

78 while (coord != NULL) { 

79 fprintf (fptr, " [%.3 f ,% . 1 5 f ] " , coord->temp , coord->p) ; 

80 if (/*m> J*/coord->next != NULL) { 

81 fprintf ( fptr , " ," ) ; 

82 } 

83 fprintf (fptr , "\n" ) ; 

84 n^ar r ay _len — 1— i ; 

85 coord= coord— >next ; 

86 } 

87 fprintf (fptr , " ] ; \ n" ) ; 

88 rr^ array _len ; 

89 if ( action . type = 'T' || action. type = 'U') { 

90 if (action. type = 'T') { 

91 fprintf (fptr, "trip:= [\n"); 

92 } else { 

93 fprintf (fptr, "q:=[\n"); 

94 } 

95 coord= holl ; 

96 while (coord != NULL) { 

97 fprintf (fptr, " [% .3 f ,% . 12 f ] " , coord->temp , coord->q) ; 

98 if (/*TO> i*/coord->next != NULL) { 
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99 fprintf (fptr , " ," ) ; 

100 } 

101 fprintf ( fptr , "\n" ) ; 

102 m=array _len — 1— i ; 

103 coord= coord— >next ; 

104 } 

105 fprintf ( fptr , " ] ; \ n" ) ; 

106 111= array_len ; 

107 if (action. type = 'U') { 

108 fprintf (fptr, "r:=[\n"); 

109 coord= hoU ; 

110 while (coord != NULL) { 

111 fprintf (fptr, " [% . 3 f ,% . 1 2 f ] " , coord->temp , coord->r ) ; 

112 if (/*m> 1*/ coord->next != NULL) { 

113 fprintf (fptr , " ," ) ; 

114 } 

115 fprintf (fptr , "\n") ; 

116 rc^ar r ay _lc n — 1— i ; 

117 coord= coord— >next ; 

118 } 

119 fprintf (fptr , " ] ; \ n" ) ; 

120 ni= array_Ien ; 

121 fprintf (fptr, " t r ipu : = [\ n" ) ; 

122 coord= hoU ; 

123 while(coord != NULL) { 

124 fprintf (fptr, " [% . 3 f , % . 1 2 f ] " , coord->temp , coord->t r ipu 

); 

125 if (/*m> _/*/coord->next != NULL) { 

126 fprintf (fptr , " ," ) ; 

127 } 

128 fprintf (fptr , "\n") ; 

129 ii^ar r ay _lc n — 1— i ; 

130 coord= coord— >next ; 

131 } 

132 fprintf (fptr , " ] ; \ n" ) ; 

133 rr^ array_len ; 

134 fprintf (fptr, " t r ipd : = [\ n" ) ; 

135 coord= holl ; 

136 while (coord != NULL) { 

137 fprintf (fptr, " [% . 3 f , % . 1 2 f ] " , coord->temp , coord->t r ipd 

); 

138 if {/*m> i*/coord->next != NULL) { 

139 fprintf (fptr , " ," ) ; 

140 } 

141 fprintf (fptr , "\n") ; 

142 m=array _len — 1— i ; 

143 coord= coord— >next ; 

144 } 

145 fprintf (fptr , " ] ; \ n" ) ; 

146 } 
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147 } 

148 fprintf (fptr, "J:=[\n"); 

149 fprintf (fptr, " %. 12 f ,% . 1 2 f ,% . 1 2 f ,%. 1 2 f " , act ion . jh , action. jv, action. jh2, 

action . j v2 ) ; 

150 fprintf (fptr , " ] ; \ n" ) ; 

151 fprintf (fptr, "Jay:= %.12f;", action. jdu); 

152 fprintf (fptr, "Jdash:= %.12f;", action. jdd); 

153 fprintf (fptr, "TripJay:= %.12f;", action. ku); 

154 fprintf (fptr, "TripJdash:= %.12f;", action. kd); 

155 fclose ( fptr ) ; 

156 } 
157 

158 void set_temp (/*/(oat temp [4] [ array -len ] , int point,*/ float min , float max) { 

159 /* A little function , to calculate the temperature points over the 

160 range between the maximum and minimum temperature levels. The random 

161 number is between and 1000, and then divided by 1000 and multiplied 

162 by the range so that we can get 200 unique points */ 

163 int i ; 

164 float range , test ; 

165 COORD * coord ; 

166 range= max— min ; //working out the range of the temperatures 

167 test= 0; //so our loop runs at least once. 

168 //To avoid divide by zero errors 

169 while (test = 0) { 

170 test= ((( float ) (rand %1000)/1000) * range)+min; 

171 } 

172 coord= hoU ; 

173 //To guarentee unique temperatures we check the current value against all 

previous values 

174 while (coord != NULL) { 

175 while (test = || test ^= coord— >temp) { 

176 tcst= ((( float ) (rand() %1000)/1000) * range)+min; 

177 coord= holl ; 

178 } 

179 coord= coord— >next ; 

180 } 

181 coord= safe_maIIoc ( sizeof (OOORD) ," Setting a temperature point"); 

182 coord— >temp= test ; 

183 coord— >next= holl ; 

184 holl= coord; 

185 } 
186 

187 void free_cond (void) { 

188 /* This frees up the memory taken by the linked list of 

189 conditions by deleting each list item one after the other 

190 from the b egining to the end. */ 

191 OOORD *del_ptr ; 

192 while (holl != NULL) { 

193 del_ptr= holl ; 

194 holl= holl->next; 
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195 free ( del_ptr ) ; 

196 } 

197 } 
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B.2.4 mach2dt.c 

1 %auto— ignore 

2 /* We link our files with a standard header */ 

3 #include " i s 1 i ng 2 d t . h" 
4 

5 void metropolis ( int chain [ length ][ height ] , float temp, ACT var , double split [3], 



6 int lattice [3] ) { 

7 /* This is the same as the previous function but now for the Union Jack 

lattice */ 

8 int i= rand { )%lcngth ; 

9 int j= rand ()%height ; 

10 int h, k; 

11 double changef= 0; 

12 double changeo= 0; 

13 double tripo ; 

14 double tripf ; 

15 double tripod : 

16 double tripfd ; 

17 double prob ; 

18 double condition ; 
19 

20 

21 changeo= energyb ( chain , var, chain [ i ][ j ] , i, j); 

22 change f= energyb (chain , var , — chain[i][j], i, j); 

23 if (var . type = 'T' ) { 

24 tripo= chain [ i ][ j ] * (( chain [( i +l)%length ][( j +l)%height ]*( chain [ i ][( j 

+ l)%height] 

25 + chain [( i+l)%length] [j ]) ) + ( chain [( i+length -l)%length ][( j+height -1)% 

height ] * 

26 (chain [ i ] [(j+height -l)%height] + chain [( i+length -l)%length ] [ j ] ) )+ (chain 

[(i+l)%length][j] 

27 * chain [ i ] [ ( j+height -l)%height ] ) +(chain [ i ] [ ( j +l)%height ]* chain [( i+ 

length-l)%length][j])); 

28 tripf= (-chain [ i ][ j ] ) * (( chain [( i +l)%length ][( j +l)%height ]*( chain [ i 

][( j+l)%height] 

29 + chain [( i+l)%length ][ j ]) ) + ( chain [( i+length -l)%length ][( j+height -1)% 

height ] * 

30 (chain [ i ] [ ( j+height -l)%height] + chain [( i+length -l)%length ] [ j ] ) )+ (chain 

[( i+l)%length] [j] 

31 * chain [ i ] [ ( j+height -l)%height ] ) +(chain [ i ] [ ( j +l)%height ]* chain [( i+ 

length-l)%length][j])); 

32 }else if(var.type = 'U' &fe ( i+j )%2 = 0) { 

33 tripo= chain [ i ] [ j ]*(( chain [( i+l)%length ] [( j +l)%height ]*( chain [( i+l)% 

length][j] 

34 + chain [ i ][( j+l)%height ]) )+ ( chain [( i+length -l)%length ][( j+height -1)% 

height ] 

35 *( chain [ i ] [(j+height -l)%height ] + chain [( i+length -l)%length ] [ j ] ) ) ) ; 

36 tripod= chain [ i ] [ j ]*(( chain [( i+length -l)%length ] [ ( j +l)%height ] 

37 *( chain [( i+length -l)%length ] [ j] + chain [ i ] [ ( j +l)%height ] ) ) 
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38 + (chain [( i+l)%length ] [ ( j+height -l)%height ] * ( chain [ i ] [( j+height -1)% 

height] + chain [( i+l)%length] [ j ]) )) ; 

39 tripf= (-chain [ i ] [ j ] ) *(( chain [( i+l)%length ] [( j +l)%height ]*( chain [( i+1) 

%length][j] 

40 + chain [ i ][( j+l)%height ]) )+ ( chain [( i+length -l)%length ][( j+height -1)% 

height ] 

41 *( chain [ i ][ ( j+height -l)%hcight ] + chain [ ( i+length -l)%lcngth ] [ j ] ) ) ) ; 

42 tripfd= (-chain [ i ] [ j ] ) *(( chain [( i+length -l)%length ] [ ( j +l)%height ] 

43 *( chain [( i+length -l)%length ] [ j] + chain [ i ] [ ( j+l)%height ] ) ) 

44 + (chain [( i+l)%length ] [ ( j+height -l)%height ]*( chain [ i ] [(j+height -1)% 

height ] + chain [( i+l)%length ] [ j ] ) ) ) ; 

45 } else if (var.type = 'U' && (i+j)%2 != 0) { 

46 tripo=chain[i][j]*((chain[i][(j + l)%height ]*chain [( i+length — l)%length 

] [ j ] ) +(chain [( i+l)%length ] [ j ] * chain [ i ] [ ( j+height -l)%height ] ) ) ; 

47 tripod= chain[i][j]* (( chain [ i ][( j +l)%hcight ]* chain [( i +l)%length ][ j ] ) 

+(chain [( i+length -l)%length ][ j ]* chain [ i ][ ( j+height -l)%height ] ) ) ; 

48 tripf= (-chain [ i ] [ j ] ) *(( chain [ i ] [( j +l)%height ]* chain [( i+length -1)% 

length ] [ j ] )+(chain [( i+l)%length ] [ j ] * chain [ i ] [ ( j+height -l)%height 

])); 

49 tripfd= (-chain [ i ][ j ]) * (( chain [ i ][( j +l)%height ]* chain [( i +l)%length ][ j 

]) +(chain [( i+length -l)%length ] [ j ]* chain [ i ] [(j+height -l)%height ] ) ) ; 

50 

51 } 

52 if(0 >= — changeo + changef) { 

53 lattice[0]= lattice[0] - 2* chain [ i ] [ j ] ; 

54 if((i+j)%2=0 &fe var.type = 'U') { 

55 lattice [1]= lattice[l] — 2* chain [ i ] [ j ] ; 

56 } else if ((i+j)%2 ! = && var.type = 'U') { 

57 lattice[2l= lattice[2] - 2* chain [ i ] [ j ] ; 

58 } 

59 if (var.type = 'T') { 

60 latticc[l]= latticc[l] — tripo + tripf; 

61 } else if (var.type ^= 'U') { 

62 lattice[3]= lattice[3]— tripo+tripf; 

63 lattice [4]= lattice [4] — tripod + tripfd ; 

64 } 

65 chain[i][j]= — chain [ i ][ j ] ; 

66 } else { 

67 prob= cxp( — (double) (l/(BOLT * temp) )*( — changeo + changef)); 

68 condition= CDND*prob ; 

69 if (rand ()903ND < condition) { 

70 lattice[0]= lattice[0] - 2* chain [ i ] [ j ] ; 

71 if((i+j)%2 = &fe var.type = 'U') { 

72 lattice[l]= lattice[l] - 2* chain [ i ] [ j ] ; 

73 } else if ((i+j)%2 !=0 &fe var.type = 'U') { 

74 lattice[2]= lattice[2] - 2* chain [ i ] [ j ] ; 

75 } 

76 if (var.type = 'T') { 

77 lattice[l]=lattice[l]— tripo + tripf; 

78 } else if(var.type = 'U') { 
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79 lattice[3]= lattice [3] — tripo + tripf; 

80 lattice[4]= lattice [4] — tripod + tripfd ; 

81 } 

82 chain [ i ] [ j ]= —chain [ i ] [ j ] ; 

83 } 

84 } 

85 split [0]= split [0] + (double) lattice [0] /( length*height ) ; 

86 if ( var . type = 'T' ) { 

87 split[l]= split [1] + (double) lattice [l]/(2*(((length-l)*(height-l))+(( 

length -l)+(height -1)) ) ) ; 

88 } else if(var.type = 'U') { 

89 split[l]= split [1] + (double) lattice [l]/(0. 5* length*height) ; 

90 split[2]= split [2] + (double) lat t ice [2] /( 0. 5 * length* height ) ; 

91 split [3]= split [3] + (double) lattice [3]/(( ( length -l)*(height -l))+(( 

length-l) + (height-l))) ; 

92 split[4]= split [4] + ( double ) 1 at t i c e [ 4 ]/((( length -1) *( height -1) ) +(( 

length -l)+(height -1)) ) ; 

93 } 

94 } 
95 

96 double energy (int chain [ length ][ height ] , ACT cond) { 

97 /* This is the energy calculation for the triangular and square lattices 

98 and returns the resultant value back */ 

99 int i , j , k , m; 

100 double total; 

101 double holdu= cond . ku ; 

102 double holdd= cond . kd ; 

103 double hu= cond.jdu; 

104 double hd= cond.jdd; 

105 double hv= cond.jv; 

106 double hh= cond.jh; 

107 double suml= 0; 

108 double sum2= 0; 

109 double sum3= 0; 

110 k= 0; 

111 for (i=0; i < length; i++) { 

112 if (k = 1) m= 3; 

113 if (k = 0) m= 0; 

114 for (j=0; j < height; j++) { 

115 if (cond. type = 'U') { 

116 switch (k) { 

117 case 0: 

118 holdd= cond.kd; 

119 holdu= 0; 

120 hd= cond.jdd; 

121 hu= 0; 

122 hh= cond . jh2 ; 

123 hv= cond. jv2 ; 

124 k= 1; 

125 break ; 
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126 case 1: 

127 holdd= 0; 

128 holdu= cond . ku ; 

129 hd= 0; 

130 hu= cond.jdu; 

131 hh= cond.jh; 

132 hv= cond.jv; 

133 k= 0; 

134 break ; 

135 default: // well you never know it could 

happen ! 

136 printf("An error has occured in the 

caIculation\n" ) ; 

137 //return (— 1) ; 

138 } 

139 } 

140 if(cond.type = 'H') { 

141 switch(k) { 

142 case 0: 

143 hv= cond . j v ; 

144 hu= 0; 

145 hd= 0; 

146 hh= 0; 

147 k= 1; 

148 break ; 

149 case 1: 

150 hv= cond . j v ; 

151 hu= 0; 

152 hd= 0; 

153 hh= cond . jh ; 

154 k= 0; 

155 break ; 

156 default: // well you never know it could 

happen ! 

157 printf("An error has occured in the 

calcuIation\n" ) ; 

158 return(-l); 

159 } 

160 } 

161 if (cond. type = 'D') { 

162 holdu= 0; 

163 holdd= 0; 

164 switch (m) { 

165 case 0: 

166 hv= cond . kd ; 

167 hu= cond . j v ; 

168 hd= 0; 

169 hh= cond . ku ; 

170 m= 1; 

171 break ; 
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172 case 1 : 

173 hv= 0; 

174 hu= cond . kd ; 

175 hd= 0; 

176 hh= 0; 

177 m= 0; 

178 break ; 

179 case 2: 

180 hv= cond . j v ; 

181 hu= 0; 

182 hd= 0; 

183 hh= cond . jh ; 

184 m= 0; 

185 break ; 

186 case 3: 

187 hv= cond .kd; 

188 hu= 0; 

189 hd= cond . j v ; 

190 hh= 0; 

191 m= 4; 

192 break ; 

193 case 4: 

194 hv= 0; 

195 hu= 0; 

196 hd= cond . kd ; 

197 hh= cond . jh ; 

198 m= 5; 

199 break ; 

200 case 5: 

201 hv= cond . j v ; 

202 hu= 0; 

203 hd= 0; 

204 hh= cond.ku; 

205 m= 3; 

206 break ; 

207 default: // well you never know it could 

happen ! 

208 printf("An error has occured in the 

calculation\n" ) ; 

209 return(-l); 

210 } 

211 } 

212 // Calculation for nearest neighbours 

213 suml= suml + ( chain [ i ][ j ]*(( hh* chain [( i +l)%length ][ j ] ) 

214 + (hv*chain [ i ] [( j+l)%height ]) + (hu* chain [( i +l)%length ][( j +1)% 

height]))) 

215 + (hd*chain [ i 1 [( j+l)%height ]*chain [( i+l)%length ] [ j ]) ; 

216 // Calculation for triplet neighbours 

217 sum3= sum3 + ( holdu* chain [ i ][ j ]*( chain [( i +l)%length ][( j +1)% 

height ] 
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218 * (chain [( i+l)%length ][ j] + chain [ i ][( j+l)%height ])) ) 

219 + (holdd*chain [ i ][( j+l)%height ]*( chain [( i+l)%length ] [j ] 

220 * (chain [( i+l)%lcngth ][( j+l)%height] + chain [ i ] [j ]))) ;; 

221 // Calculation for sites 

222 sum2= sum2 + chain [ i ] [ j ] ; 

223 } 

224 if(i%2 = 0) { // to allow the intial diagonal of the row to alternate 

225 k= 1; 

226 } else { 

227 k= 0; 

228 } 

229 } 

230 // adding the three sums together and applying required factors 

231 totaI= (— sum3)— suml — ( cond . h*sum2) ; 

232 // returning the result 

233 return total; 



234 } 

235 
236 

237 double energyb(int chain [ length ][ height ] , ACT cond, int change, int i, int j) { 



238 /* This is the energy calculation for the triangular and square lattices 

239 and returns the resultant value back */ 

240 int k, m; 

241 double total ; 

242 double holdu= cond . ku ; 

243 double holdd= cond . kd ; 

244 double holdua= holdu ; 

245 double holdda= holdd ; 

246 double hu= cond.jdu; 

247 double hd= cond.jdd; 

248 double hv= cond.jv; 

249 double hh= cond.jh; 

250 double hua= hu ; 

251 double hva^ cond.jv2; 

252 double hda^ hd ; 

253 double hha= cond.jh2; 

254 double suml= 0; 

255 double sum2= 0; 

256 double sum3= ; 

257 if (( i+j )%2 = 0) { 

258 k= 1; 

259 } else { 

260 k= 0; 

261 } 

262 if(i%2 = &fe cond . type = 'D') { 

263 switch (j%3) { 

264 case 0: 

265 0; 

266 break ; 

267 case 1: 
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268 m= 1; 

269 break ; 

270 case 2: 

271 m= 2; 

272 break ; 

273 default : 

274 printf{"An error has occured"); 

275 exit(-l); 

276 } 

277 } else if(cond.type = 'D') { 

278 switch (j%3) { 

279 case 0: 

280 m=3; 

281 break ; 

282 case 1: 

283 m= 4; 

284 break ; 

285 case 2: 

286 m= 5; 

287 break ; 

288 default : 

289 printf("An error has occured"); 

290 exit(-l); 

291 } 

292 } 
293 

294 if(cond.type = 'U') { 

295 switch (k) { 

296 case 0: 

297 holdd= 0; 

298 holdda= cond . kd ; 

299 holdu= 0; 

300 holdua= cond . ku ; 

301 hd= 0; 

302 hda^ 0; 

303 hu= 0; 

304 hua= 0; 

305 hh= cond . jh2 ; 

306 liv= cond . j v2 ; 

307 hha= cond . jh ; 

308 hva^ cond.jv; 

309 k= 1; 

310 break ; 

311 case 1: 

312 holdd= cond.kd; 

313 holdda= 0; 

314 holdu= cond . ku ; 

315 holdua= 0; 

316 hd= cond.jdd; 

317 hda^ cond.jdd; 
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318 hu= cond.jdu; 

319 hua^ cond.jdu; 

320 hh= cond.jh; 

321 hv= cond . j v ; 

322 hha^ cond. jh2 ; 

323 hva^ cond.jv2; 

324 k= 0; 

325 break ; 

326 default: // well you never know it could happen! 

327 printf("An error has occured in the calculation \n" ) ; 

328 //return(-l) ; 

329 } 

330 } 

331 if (cond. type = 'H') { 

332 holdd= 0; 

333 holdda= 0; 

334 holdu= 0; 

335 holdua= 0; 

336 switch(k) { 

337 case 0: 

338 hu= 0; 

339 hua^ 0; 

340 hd= 0; 

341 hda^ 0; 

342 hh= 0; 

343 hha= cond.jh; 

344 k= 1; 

345 break ; 

346 case 1: 

347 hu= 0; 

348 hua^ 0; 

349 hd= 0; 

350 hda= 0; 

351 hh= cond . jh ; 

352 hha^ 0; 

353 k= 0; 

354 break ; 

355 default: // well you never know it could happen! 

356 printf("An error has occured in the calculation \n" ) ; 

357 return(-l); 

358 } 

359 } 

360 if (cond. type = 'D') { 

361 holdu= 0; 

362 holdua= 0; 

363 holdd= 0; 

364 holdda= 0; 

365 switch (m) { 

366 case : 

367 hv= cond.kd; 
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368 hva^ cond . kd ; 

369 hu= cond.jv; 

370 hua= cond.jv; 

371 hd= 0; 

372 hh= cond.ku; 

373 hha^ cond . jh ; 

374 m= 1; 

375 break ; 

376 case 1: 

377 hv= 0; 

378 hu= cond . kd ; 

379 hua^ cond . kd ; 

380 hd= 0; 

381 hh= 0; 

382 hha= cond . ku ; 

383 m= 0; 

384 break ; 

385 case 2: 

386 hv= cond . j v ; 

387 hva^ cond.jv; 

388 hu= 0; 

389 hd= 0; 

390 hh= cond.jh; 

391 0; 

392 break ; 

393 case 3: 

394 hv= cond.kd; 

395 hva= cond . kd ; 

396 hu= 0; 

397 hd= 0; 

398 hh= 0; 

399 hha^ cond.ku; 

400 m= 4; 

401 break ; 

402 case 4: 

403 hv= 0; 

404 hu= 0; 

405 hd= cond.kd; 

406 hda= cond . kd ; 

407 hh= cond . jh ; 

408 5; 

409 break ; 

410 case 5: 

411 hv= cond.jv; 

412 hva^ cond.jv; 

413 hu= 0; 

414 hd= cond . kd ; 

415 hda^ cond . kd ; 

416 hh= cond . ku ; 

417 hha^ cond.jh; 



B.2 Code listings for main program 



131 



418 m= 3; 

419 break; 

420 default: // well you never know it could happen! 

421 printf("An error has occured in the calculation \n" ) ; 

422 return (-1); 

423 } 

424 } 

425 // Calculation for nearest neighbours 

426 suml= suml + ( change *(( hh* ( chain [( i +l)%length ][ j ] ) 

427 + (hv*chain [ i ] [ ( j+l)%hcight ] ) + ( hu* chain [( i +l)%length ][( j +l)%height ] ) 

428 + (hd*chain [( i+l)%lcngth ] [ ( j+(hcight -l))%height ] ) 

429 + (hha*chain [( i+(length -l))%length ] [ j 1 ) + (hva* chain [ i ][( j +(height -1) )%height 

]) 

430 + (hua*chain [( i+(length -l))%length ] [ ( j+height -l)%hcight ] ) 

431 + (hda*chain [( i+(length -l))%length ] [( j +l)%height ] ) ) ) ) ; 

432 // Calculation for triplet neighbours 

433 sum3= sum3 + ( change *(( holdu* chain [( i +l)%length ][( j +l)%height ] 

434 * (chain [( i+l)%length ][ j] + chain [ i ][( j+l)%height ]) ) 

435 + (holdu*chain [( i+length -l)%length ][( j+height -l)%height ] 

436 * (chain [( i+length -l)%length ] [ j] + chain [ i ] [(j+height -l)%height ]) ) 

437 + (holdd*chain [( i+l)%lcngth ][( j+height -l)%hcight ] 

438 * (chain [ i ] [( j+hcight - l)%hcight ] + chain [( i+l)%length ] [ j ] ) ) 

439 + (holdd*chain [( i+length -l)%length ][( j +l)%height ] 

440 * (chain [ i ][( j+l)%height] + chain [( i+length -l)%length ][ j ]) ) 

441 + (holdda*chain [( i+l)%lcngth ][ j ]* chain [ i ][( j +l)%height ] ) 

442 + (holdua* chain [( i+l)%length ] [ j ]* chain [ i ] [(j+height -l)%height ] ) 

443 + (holdda* chain [( i+length -l)%length ] [ j ]* chain [ i ] [ ( j+height -l)%height ] ) 

444 + ( holdua* chain [ ( i+length -l)%length ] [ j ] * chain [ i ] [ ( j +l)%height ] ) ) ) ; 

445 // Calculation for sites 

446 sum2= sum2 + change ; 

447 // adding the three sums together and applying required factors 

448 total= ( — sum3)— suml — ( cond . h*sum2) ; 

449 // returning the result 

450 return total; 



451 } 
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Programming II 
Mean Field Approximation 



C.l Programmers guide 

The broad purpose of this code is to perform mean field simulations of the Ising model 
on various lattices. As with the last program, it does this by asking the user initially 
for the conditions for the system that is to be studied. These include the shape of the 
lattice, the temperature range to be studied over and the interaction strengths of the 
lattice. It then performs simulations of systems at temperatures at regular intervals 
along the temperature range. The program then outputs the results to a file which can 
be read into a mathematical program. 

C.1.1 Flow 

Initially the user is asked for the initial conditions of the system they wish to study, 
and then the temperature range they would like to study. After these have been 
stored, the user is then asked how many data points they require. After this the 
program calculates each data point by iteratively using the formulas from Chapter [5] 
until either the value of the sublattice magnetisation is constant, or 10,000 iterations 
have taken place, whichever occurs first. The value of magnetisation is then recorded for 
both the partially uncoupled and coupled equations along with the temperature. The 
program then moves on to the next data point, using the previous magnetisation value 
as the initial value. The temperature values are calculated at regular intervals in the 
temperature range. After the required number of data points have been calculated, the 
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results are output to a data file which can then be read into a mathematical program. 
Again the simulation time is given on the screen after the program has finished. 



C.2 Code listings 
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Code listings 



C.2.1 



meanl.h 



1 %auto— ignore 

2 ^include <stdio.h> 

3 #include <stdlib.h> 

4 #include <math . h> 

5 #include <time.h> 
6 

7 #define BOLT 8.617343*pow(10 , -5) /*The boltzman contant*/ 
8 

9 typedef struct coordinate { 

10 long double sigma ; 

11 long double mu; 

12 long double unsigma ; 

13 long double unmu; 

14 float temp; 

15 struct coordinate *next; 

16 } (X)ORD; 
17 

18 double set up _gcn ( char *, char *); 

19 void * saf e _malloc ( s iz e _t , char *); 

20 void *safe_realIoc (int *, size_t , char *) ; /* Does the same for realloc */ 

21 void write_2d_array (char *, double, double, double, char); 

22 void frcc_cond (void); /* Frees the memory taken by linked list */ 

23 long double ltanh(long double); 
24 

25 extern OOORD *holl ; 
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C.2.2 meanl.c 

1 %auto— ignore 

2 #include "meanl.h" 
3 

4 OOORD *holl= NULL; 
5 

6 int main() { 

7 int i , timeleft , points , seconds , minutes , hours , c , d; 

8 time_t timer ; 

9 OOOm * coord ; 

10 float min , max, minmax, percent , interval , T, q; 

11 char type ; 

12 double J, K, B; 

13 long double x, y, u, v, beta, m, n, a, b; 

14 char filename [ 1 0] ; 
15 

16 x= 0; 

17 y= 0; 

18 while(type != 'S' &fe type != 'T' &fe type != 'U' &fe type != 'N' &fe type != 'H' 

&&: type != 'D') { 

19 printf ("What type of lattice would you like to study? \n ([S]quare, [T 

] riangular , [U] nion Jack , [N] ext Nearest , [H] exagonal , [D] ouble 
Hexagonal ) " ) ; 

20 fscanf ( stdin , "%ls" , &type); 

21 type= toupper { type ) ; 

22 } 

23 switch (type) { 

24 case ' S ' : 

25 q= 2; 

26 break ; 

27 case 'T': 

28 q= 3; 

29 break ; 

30 case 'U': 

31 q= 2; 

32 break ; 

33 case 'H': 

34 q = l-5; 

35 break ; 

36 default : 

37 p ri n t f (" There has been an error"); 

38 exit(-l); 

39 } 

40 if(type = 'U') { 

41 J= setup_gen (" interaction of the square lattice", "real"); 

42 K= setup_gen (" interaction of the diagonal lattice", "real"); 

43 } else { 

44 J= setup_gen (" interaction of nearest neighbours", "real"); 

45 } 
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46 B= setup_gen (" strength of overal field", "real"); 

47 x= (long double) setup_gen (" intial value for sigma" , "real, less than or equal 

to +/- 1") ; 

48 if(type = 'U') { 

49 y= (long double) setup_gen (" intial value for mu" , "real, less than or 

equal to +/- 1" ) ; 

50 } 

51 min= ; 

52 max= min ; 

53 while (min = max) { 

54 min= ( float ) set up _gen (" minimum temperature", "real above zero"); 

55 ma;^ ( float ) setup_gen ("maximum temperature", "real above zero"); 

56 if(min < 0) min= 0; //As we are working in Kelvin it makes little 

sense to have negative numbers 

57 if (max < 0) max= 0; 

58 } 

59 //to save retyping if the temperatures are inputted the wrong way round 

60 if (min > max) { 

61 minmax= min ; 

62 min= max; 

63 max= minmax ; 

64 } 

65 points= ( int ) setup-gen (" data points", "integer above zero"); 

66 p r i nt f (" Enter output filename: "); 

67 if (fscanf (stdin , "%s" , ^filename ) != NULL) { 

68 printf ("\n") ; 

69 fflush (stdin) ; 

70 } 

71 timer= time (NULL) ; // starting the timer 

72 // creating our temperature array 

73 interval= (max— min) / points ; 

74 c= 0; 

75 d= 0; 

76 coord= holl ; 

77 coord= safc_malloc ( sizeof (OCXDRD) , "New data point"); 

78 coord— >temp= min; 

79 coord— >sigma= x; 

80 coord— >mv^ y; 

81 coord— >unsigma= x; 

82 coord— >immu= y; 

83 coord— >next= holl ; 

84 holl= coord; 

85 u= x ; 

86 v= y; 

87 m= 10; 

88 n= 10; 

89 a^ 10; 

90 b= 10; 

91 for(i= 1; i < points; i++) { 

92 T= min+( i * in t er V al ) ; 
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93 beta= (long double )( 1 / (BOLT * T) ) ; 

94 while (x != m | | y != n) { 

95 m= X ; 

96 n= y; 

97 if(type = 'U') { 

98 x= beta*((2*K*m)+(2*J*n)+B) ; 

99 y= beta*((2* J*m)+B) ; 

100 } else { 

101 x= beta * ( ( q* J*m)+B) ; 

102 } 

103 if(c = 0) { 

104 x= Itanh (x) ; 

105 y= Itanh (y) ; 

106 } else { 

107 x= (Itanh (x)+ m) /2; 

108 y= (ltanh(y)+ n)/2; 

109 } 

110 C++; 

111 if (c = pow(10 ,8) ) break; 

112 } 

113 c= 0; 

114 while(u != a || v != b &fe type = 'U') { 

115 a^ u ; 

116 b= v; 

117 if(d= 0) { 

118 u= ltanh(beta*((2*K*a)+B)) ; 

119 v= ltanh(beta*(4*J*a)+B) ; 

120 } else { 

121 u= (Itanh (beta*((2*K*a)+B) )+a) /2; 

122 v= (Itanh (beta*((4*J*a)+B) )+b) /2; 

123 } 

124 d++; 

125 if (d = pow(10 ,8) ) break; 

126 } 

127 d= 0; 

128 coord= safe_malloc ( sizeof (OOORD) , "New data point"); 

129 coord ->temp= T; 

130 coord— >sigma= x; 

131 coord— >mi:^ y; 

132 coord— >unsigma= u; 

133 coord— >uninu= v; 

134 coord— >next= holl ; 

135 holl= coord ; 

136 m= 10; 

137 n= 10; 

138 a^ 10; 

139 b= 10; 

140 percent= (float) ( i+l)/points ; 

141 seconds= ( time (NULL) — timer ) ; 

142 timeleft= (long int )((( time (NULL)— timer )/ percent ) — seconds ) ; 
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143 seconds= timeleft %60; 

144 minutcs= abs ( t i m c 1 c f t /60) %60; 

145 hours= abs(timeleft/(60*60)) ; 

146 printf ("%.2f %% done , time left: %.2d hours, %.2d minutes, %.2d seconds 

\r" , percent *100 , hours, minutes , seconds ) ; // To show the program is 
still working 

147 } 

148 write_2d_array ( filename , J, K, B, type); 

149 frcc_cond () ; 

150 printf ("\n Now run maple to plot graph\n" ) ; // handy hint to the user 

151 seconds= ( time (NULL)-timer ) %60; 

152 minutes= ( abs (( time (NULL)— timer— seconds ) /60) ) %60; 

153 hours= abs { ( time (NULL)— timer — seconds — (minutes*60)) / (60*60)) ; 

154 print f (" Runtime : %ld hours, %ld minutes, %ld seconds\n" , hours, minutes, 

seconds); // how long it too 

155 return 0; 

156 } 
157 

158 long double ltanh(long double x) { 

159 if (fabsl(x) < powl (2.0 , ( - 55 .0) ) ) { 

160 x= x; 

161 } else { 

162 x= tanhl (x) ; 

163 } 

164 return x ; 

165 } 
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C.2.3 setup. c 

1 %auto— ignore 

2 /* We link our files with a standard header */ 

3 #include " isling2 dt . h" 
4 

5 void * safe_malloc (size_t size , char *location) { 

6 /* This function performs a malloc operation but checks to 

7 see if there is enough free memory for this to be written to. */ 

8 void * ptr ; 

9 ptr= malloc ( size ) ; 

10 if (ptr = NULL) { 

11 printf("Out of memory at function: %s\n" , location ) ; 

12 exit(-l); 

13 } 

14 return ptr ; 

15 } 
16 

17 void * safe_realloc ( int *array , size_t size, char *location) { 

18 /* This is like the safc-malloc but for realloc. */ 

19 void * ptr ; 

20 ptr= realloc ( array , size); 

21 if (ptr = NULL) { 

22 printf("Out of memory at function: %s\n" , location ) ; 

23 exit(-l); 

24 } 

25 return ptr ; 

26 } 
27 

28 int setup_chain (char *quant) { 

29 /* This function is a left over from a previous program for 

30 the one dimensional Ising model. A two dimensional version 

31 is being looked into.*/ 

32 int num; 

33 printf {" Enter the %s of chain required: (0,1,...) ", quant); 

34 fscanf ( stdin , "%d" , &num) ; 

35 printf ("\n") ; 

36 return num; 

37 } 
38 

39 void i ni t i alise _c hai n ( int chain [ length ][ height ] ) { 

40 /* This function take the chain and set all site in the "UP" 

41 position */ 

42 int i , j ; 

43 for(j=0; j < height; j++){ 

44 for(i=0; i < length; i++){ 

45 if(chain[i][j] !=UP) chain [ i ][ j ]= UP; 

46 } 

47 } 

48 } 
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49 

50 double setup_gcn (char *item, char * units) { 



51 /* This is the basic input function , as so many of our variables 

52 have similar questions I decided to write a function for them*/ 

53 double i ; 

54 p r i nt f {" Enter the value of the %s : (%s) ", item , units ) ; 

55 if (fscanf (stdin , "%lf", &i ) ) { 

56 printf ("\n") ; 

57 } else { 

58 printf {" \n EFiROR: No input found. \n" ) ; 

59 clearerr(stdin); 

60 fflush (stdin) ; 

61 setup_gcn ( item , units ) ; 

62 fflush (stdin) ; 

63 } 

64 return i ; 



65 } 
66 

67 void write_2d-array (char output [] , ACT action) 

68 /* Take the arrays in x arid y both of length array-ten and write a file in 

69 maple output format to the file named filename */ 

70 { 



71 int i ; 

72 int rr^ array_len ; 

73 OOOFiD * coord; 

74 FILE *fptr ; 

75 coord= holl ; 

76 fptr= fopen (output, "w" ) ; 

77 fprintf (fptr , "p:=[\n" ) ; 

78 while (coord != NULL) { 

79 fprintf (fptr, " [% . 3 f , % . 1 5 f ] " , coord->temp , coord->p) ; 

80 if (/*m> J*/coord->next != NULL) { 

81 fprintf (fptr , " ,") ; 

82 } 

83 fprintf ( fptr , "\n" ) ; 

84 rr^array -len — 1— i ; 

85 coord= coord— >next ; 

86 } 

87 fprintf (fptr , " ] ; \ n" ) ; 

88 array_len; 

89 if ( action . type ^= 'T' || action. type = 'U') { 

90 if (action . type = 'T') { 

91 fprintf (fptr, "trip:= [\n"); 

92 } else { 

93 fprintf (fptr, "q:=[\n"); 

94 } 

95 coord= holl ; 

96 while (coord != NULL) { 

97 fprintf (fptr, " [% . 3 f , % . 1 2 f ] " , coord->temp , coord->q) ; 

98 if i/*m> i*/coord->next != NULL) { 
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99 fprintf (fptr , " ," ) ; 

100 } 

101 fprintf (fptr , " \n" ) ; 

102 iii=array_len — 1— i ; 

103 coord= coord— >next ; 

104 } 

105 fprintf ( fptr , " ] ; \ n" ) ; 

106 m= array_len ; 

107 if ( action . type = 'U') { 

108 fprintf (fptr, "r:=[\n"); 

109 coord= lioll ; 

110 while (coord != NULL) { 

111 fprintf (fptr, " [% . 3 f , % . 1 2 f ] " , coord->temp , coord->r ) ; 

112 if (/*m> 1*/ coord->next != NULL) { 

113 fprintf (fptr , " ," ) ; 

114 } 

115 fprintf (fptr , "\n") ; 

116 rr^array _Ien — 1— i ; 

117 coord= coord— >next ; 

118 } 

119 fprintf (fptr , " ] ; \ n" ) ; 

120 111= array_len ; 

121 fprintf (fptr, " t ripu : = [\ n" ) ; 

122 coord= holl ; 

123 while (coord != NULL) { 

124 fprintf (fptr, " [% .3 f ,% . 1 2 f ] " , coord->temp , coord->t r ipu 

); 

125 if (/*m> J*/coord->ncxt != NULL) { 

126 fprintf (fptr , " ," ) ; 

127 } 

128 fprintf (fptr , "\n") ; 

129 rr^array -len — 1— i ; 

130 coord= coord— >next ; 

131 } 

132 fprintf ( fptr , " ] ; \ n" ) ; 

133 ni= array_Ien ; 

134 fprintf (fptr, " t r ipd : = [\ n" ) ; 

135 coord= holl ; 

136 while (coord != NULL) { 

137 fprintf (fptr, " [%.3 f ,%. 1 2 f ] " , coord->temp , coord->t r ipd 

); 

138 if (/*TO> _?*/coord->ncxt != NULL) { 

139 fprintf ( fptr , " ," ) ; 

140 } 

141 fprintf (fptr , "\n") ; 

142 iii=array_len — 1— i ; 

143 coord= coord— >next ; 

144 } 

145 fprintf (fptr , " ] ; \ n" ) ; 

146 } 



C.2 Code listings 



143 



147 } 

148 fprintf (fptr, "J:=[\n"); 

149 fprintf (fptr, " %. 12 f , % . 1 2 f , % . 1 2 f , % . 1 2 f " , act ion . jh , action. jv, action. jh2, 

action . j v2 ) ; 

150 fprintf (fptr , " ] ; \ n" ) ; 

151 fprintf (fptr, "Jay;= %.12f;", action. jdu); 

152 fprintf (fptr, "Jdash:= %.12f;", action. jdd); 

153 fprintf (fptr, "TripJay:= %.12f;", action. ku); 

154 fprintf (fptr, "TripJdash:= %.12f;", action. kd); 

155 fclose ( fptr ) ; 

156 } 
157 

158 void set -temp ( /* /Z o a f temp [ 4 ][ array _len ] , int point,*/ float min , float max) { 

159 /* A little function , to calculate the temperature points over the 

160 range between the maximum and minimum temperature levels. The random 

161 number is between and 1000, and then divided by 1000 and multiplied 

162 by the range so that we can get 200 unique points */ 

163 int i ; 

164 float range , test ; 

165 COORD * coord; 

166 range= max— min ; //working out the range of the temperatures 

167 test= 0; //so our loop runs at least once. 

168 //^o avoid divide by zero errors 

169 while (test = 0) { 

170 test= ((( float ) (rand %1000)/1000) * range)+min; 

171 } 

172 coord= hoU ; 

173 //^o guarentee unique temperatures we check the current value against all 

previous values 

174 while (coord != NULL) { 

175 while (test = || test = coord— >tcmp) { 

176 test= ((( float ) (rand %1000)/1000) * range)+min; 

177 coord= hoU ; 

178 } 

179 coord= coord— >next ; 

180 } 

181 coord= safe_malIoc ( sizeof (OOORD) ," Setting a temperature point"); 

182 coord— >tcmp= test; 

183 coord— >next= holl ; 

184 holl= coord; 

185 } 
186 

187 void free_cond (void) { 

188 /* This frees up the memory taken by the linked list of 

189 conditions by deleting each list item one after the other 

190 from the begining to the end. */ 

191 OOORD *del_ptr ; 

192 while (holl != NULL) { 

193 del_ptr= holl ; 

194 hoIl= holl->next; 
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195 free (del.ptr) ; 

196 } 

197 } 



